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ABSTRACT 


The current distribution and impedance of a thin cylindrical 
antenna with parallel orientation to the static magnetic field of a 
lossy magnetoplasma is calculated with the method of moments. The 
electric field produced by an infinitesimal current source is first 
derived. Results are presented for a wide range of plasma parameters. 
Reasonable answers are obtained for all cases except for the overdense 
hyperbolic case. A discussion of the numerical stability is included 
which not only applies to this problem but other applications of the 
method of moments . 
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CHAPTER I. INTRODUCTION 


Since man has launched space vehicles into the ionosphere, there 
has been strong interest in the effect of plasmas on antennas. A 
primary interest is understanding the effects of the plasma on the- 
current distribution and the antenna impedance. Since cylindrical 
antennas have a simple geometry, are well-understood in free space, and 
have wide usage, they are most commonly studied in plasma problems. 

Because one cannot analytically solve Maxwell's equations for a lossy 
magnetoplasma with cylindrical boundary conditions, a .current distribution 
usually is assumed. This assumed current distribution, though, represents 
only the limiting case of an infinitesimal antenna. Approximations made 
by other workers give results valid for inf inite length antennas . These 
approximations simplify the analytical work but give results of questionable 
validity for the finite antenna. 

The approach used in this report eliminates these approximations. 

The method of moments is used to solve for current distribution on the 
antenna. Basically, this method converts an integral equation into a 
matrix equation which can be easily solved by computer methods. From the 
current distribution one can trivially find the impedance at the .feed 
point. Thus, the .current distribution and impedance for the finite 
antenna over a wide range of plasma densities and' applied magnetic fields 
can be calculated. . ..... 

The geometry used in this problem is shown in Figure 1. The 
orientation of the antenna is restricted to the z direction .parallel to 
the applied magnetic field. The antenna is described by the half- 
length H and radius a. For a thin-wire antenna, the condition a << H 
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Figure 1. Geometry of Problem 
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must hold. For the method of moments, the antenna is divided into an 
odd number of segments N. 

The first step in applying the method of moments is solving for 
the electric field due to an infinitesimal source. Exact details of 
how the method of moments is applied will be discussed. Results are 
given for a wide range of plasma parameters and are compared to the 
limiting cases of free space and quasi-static models. Comparisons are 
also made to experimental data. A sensitivity theory is discussed 
which applies not only to this problem but to many other applications 


of the moment method 
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CHAPTER II. FIELDS DUE TO INFINITESIMAL SOURCE 
2.0 An Overview 

To determine the current distribution along an antenna it is first 
necessary to find the fields produced by an infinitesimal current 
source. The derivation given here used the notation and technique 
developed by Mittra and Deschamps [1], In their paper they outline the 
derivation of the field components for any orientation of an 
infinitesimal source. The parallel orientation case worked out in 
detail here, however, is not the same as the specific case completed 
in their paper. 

In this report the plasma will be characterized by normalized 


plasma parameters X, Y 

, and 

Z where 



uj 2 


.10 1 2 


X = ( 


Y 2 = 

- . 

, = v . 

1 

. CO / * 


\ to 1 

0) 

r 

N e 2 1 

1/2 



0) = • 
P 

o 

me 

= 

plasma 

frequency 


o J 


eB 

oj c = — = cyclotron frequency . 

a) = operating frequency 

y = collision frequency 

N q = average electron density 

e = magnitude of electronic charge 

m = mass of electron 

B q = applied magnetic field. 

The^use - of _ X7~Y, and - Z~ parameters allows one to indicate results on 
commonly used plasma diagramfe. 
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Before giving the lengthy details, we wish to give an overview- 
of the derivation. Maxwell's equations can be combined to yield 
a relation for the electric field caused by an infinitesimal source 
in the plasma. Taking the Fourier transform allows one to easily solve 
for the transform of the z component of the electric field. The main 
feature of this technique is that in performing the inverse Fourier 
transform, the integral can be broken into a sum of functions . such 
that one can analytically integrate terms that are singular when 
observation and source points coincide. These terms represent the 
near field and are dominant in the determination of the antenna current 
and impedance. Terms that are finite when observation and source 
points coincide need to be numerically integrated, but these represent 
far-field terms and need not be as accurate because of their small 
relative size. The rest of this chapter gives the details; the final . 
expressions for the field due to the infinitesimal source are given 
by Equations (2.80) through (2.88). 

2.1 Maxwell's Equations 

Using e'' a>t time convention. Maxwell's equations for the plasma 

are 


curl E = -jrnyH (2.2a) 

curl H = jwe o e • E + J (2.2b) 

div H = 0 (2 . 2c) 

div e • E = p/e (2. 2d) 

. o 


where 
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-je 


e = 




(2.3) 


U = 1 - jZ 
e = 1 - UX/(U 2 - Y 2 ) 


YX 


2 2 
U - Y 


e = 1 - X/U. 
z 


When (2.2a) and (2.2b) are combined, the following equation is 
obtained 

_ 2 — — — 2 2 
curl curl E = k e • E - jwyJ, k = to ye . (2.4) 

o o o 

At this point it becomes convenient to transform the problem 

into k-space by applying the Fourier transform. Let E' and J' 

represent the transforms of E and J. , The transform E'. of E is 

x x 


defined by 

00 j (k x+k y+k z ) 

E ' (k , k , k ) = / / / E„(x, y, zj e x y Z dxdydz 

A A y Z A 

J —00 

(2.5) 

with the transforms of other components defined similarly. Note that 
the del operator can be replaced by jk. Then, using a vector 
identity, the transform of (2.4) can be written as 


k(k • E ? ) - k *k E' + k 2 e E' = jtoyJ' 

o 


( 2 . 6 ) 


Next one can expand this equation by components and write it in matrix 
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form 



— 


— 


— 


— 


— 



r 

e 

-je' 

0 


2 2 
-k-k 

k k 

k k 

I 

E' 


J' 





y z 

x y 

X z 


X 


X 

l -k 2 

je' 

e 

0 


k k 

2 2 
-k-k 

k k 

/ 

^E* 

=-jojp 

J' 

^ ° 




x y 

x z 

y z 


y 

. y 


0 

0 

e 


k k 

k k 

2 2 
-k-k 


E' 


j’ 



z 


x z 

y z 

x y 


; z 


z 

i 

— 


— 


— 


— 

— 


— — 


(2.7) 


Because of the radial symmetry about the z axis, k. can be represented 
by spherical coordinates (T, 4 1 , a), 

k = T sinfcos a, k = T sin V sin a, k = V cos (2.8) 

x y z 

Making these changes and collecting terms, we obtain 


ME' = j ^ J' 


(2.9) 


M = 


e-r^(l - sin^V cos^a) 


2 2 t. t. t. 

je'+r^ sin f sin a cos a e-T^(l - sin^'l' sin^a) 


2 2 2 
-je’+r sin V sin a cos a T sin ¥ cos ^ cos a 


, 2 „ 


sin f cos f cos a 


2 

F^sin 'H cos ¥ sin a 


2 

r^sin i cos sin a 


2 2 

e - T.sin f 
z 1 


where T, = ~r/k . 

1 o 

To solve for E’ we can write 

E' = M -1 j J'. (2.10) 

o 

Finding M ^ is a straightforward but lengthy algebraic manipulation. 
One finally obtains 
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M -1 = [a i j ] /A (2.11) 

A = det(M) = rf (e sin^¥ + e cos^¥) - [e e (1 + cos^¥) + (e^ - e 1 ^) sin^¥] 
1 z 1 z 

? 2 

+ (e - e* )e 

z 

p 4 • 2 .... 2 „2 r . 2 , ■ . 2 ... . 2 . . , 

a.. = r.sin ¥ cos a - [e sin ¥ + e (1 - sm ¥ sm a)J + e e 

11 i 1 z z 

A 2 2 2 

a 01 = r, sin ¥ sin a cos a - F n sin ¥(e sin a cos a -je f ) - je’ e 

21 1 1 z J/J z 

2 2 

a^ = r sin ¥ cos ¥ [T^ cos a - (e cos a - je* sin a)] 

4 2 2 2 

a 10 = T sin V sin a cos a - ..It sin . V (e sin a cos a + je') + je' e 

12 1 1 z J/J z 

a 0 „ = rj sin^ ¥ sin^ a - T^[e sin^ ¥ + e (1 - sin^ ¥ cos^ a)] + e e 

22 1 l l z z 

2 2 

a 32 .= r x sin ¥ cos ¥ [T^ sin a - (e sin a + je’ cos a)] 


a 13 * F 1 

sin ¥ cos ¥ 

2 

[T cos a - 

(e 

cos a + je' 

sin a) ] 

a = 

23 1 

sin ¥ cos ¥ 

[T^ sin a - 

(e 

sin a - je' 

cos a)] 

a 33 = r i 

2 2 
cos ¥ - T ± 

e(l + cos^ 

¥) 

. 2 ,2 
+ e - e 



where T = r/k . 

1 o 

It should be noted that det(M) = 0 is the dispersion relation for 

the plasma. Since, for a fixed polar angle ¥, this is a second-degree 
2 

polynomial in there will be exactly two roots. These two roots are 
the squares of the indices of refraction for two characteristic plane 

waves traveling at an angle ¥ to the magnetic fields _The_p_roper _square 

root must be taken to satisfy the radiation condition,. One can appreciate 
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the conceptual .difficulties inherent in understanding plasma problems 
since unlike free space, there is not even a unique wavelength in a 
given direction. The exact form of the indices of refraction will be 
given shortly. , 

To solve our problem, only one of the nine terms of (2.11) will be 
needed. For parallel orientation of the antenna, there is only source 
current in the z direction. Only is constrained by the boundary 
condition. Therefore, the rest of this chapter is directed to finding 
the inverse transform of 

e ; ■ (a 33 /4) k 1 j ; • (2 - i2 > 

o 

2.2 Inverse Transform 


Conceptually all that remains is to perform the inverse transform 

oo 

3 


e 2 <*. y. - f- ( ) 

o 



< ) e' 2k " r dk. 


(2.13) 


Recall that is unity since our problem is formulated for a Dirac <5 
source. If spherical coordinates are also used for the observation point, 
(2.13) becomes 


00 . ,77 2 If 

E CR. e, ♦) - & ( i ) f (f (-¥ )e -jrR[ sin'SCiM 
° 0 v 0 


-cos 0 cos f] 


• r sin V da d¥ dr 


where x = R sin 0 cos V . , y = R sin. 0 sin d> , z = R cos 0.. 

First we can relatively easily perform the a integration using the 
following standard integration formula [2], 
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2tt 


J (fR sin 0 sin V ) = 




-jTR sin 8 : sin f cos(a^<j 


da. (2.15) 


0 

Applying (2.15) to (2.14) produces 


" TT 


‘P<£> J J ( “T } -o' 
Vo 


(r R sin 6 sin 40 e cos 0COS ^ r 2 sin V d'f dr 


(2.16) 


a 33 rjcos 2 f - r 2 e( 1 + cos 2 '!') + e 2 - s’ 2 

T^(e sin 2 f + e cos 2 '!) - r 2 [e e (1 + cos 2 40 + (e 2 - E ,2 sin 2 'i , ]+(e 2 -e ,2 )e 
1 z 1 z z 

r ! - r/ V 

At this point it is convenient to write A in factored form. 

A(r, 4 1 ) = (e sin 2 4 ( + e z cos 2 4')(r 2 - n 2 ) (r 2 - n 2 ) (2.17) 

where 

2 2 B(4Q ± / B 2 (4Q - 4 A(f) CQQ 

n l’ n 2 2A(4') 

A(4>) = e sin 21 ! + e cos 2 V 

z 

. ■ B(4') = e e (1 + cos 2 4') + (e 2 - e’ 2 ) sin 2 4' 

z 

C (40 = (e 2 - e' 2 ) e ■ . 

z 

Recall n^ and n ^ represent indices of refraction in the plasma mentioned 
earlier. 

The ability to do the integration: in (2.16) depends upon it being 

a 33 

possible to break into a sum of functions with different dominant 

4 

orders - of - r. AFlnentioned in the overview, it will be shown that those 
terms which are singular when observation and source point coincide can 
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be Integrated in closed form. These represent near-field terns and 
will be dominant in this problem. Those terms that are finite when 
observation and source points coincide will need to be integrated . 
numerically but will be small compared to singular terms. 

By successive synthetic division one can show 


a„ 9 _ _ _ 2 m r 2 k 2 (e 2 + £ ' 2 cos 2 4')sin 2 'F. 

33 „2 _ cos T o 

A " AOO .2 


(2.18) 




2 2 
- k sin 1 
o 


2 . 2 r , 2 , 2 . , 2 , 2 W 2 , ,2 2 w '. ^2 2 / 2 , ,2 X . 4 , 

(r k Q [e(e -e' ) - (n^+ n^) (e +e ' cos 4')] +n^ ^ (e +e cos; Tjk o > 

2 2 2 2 ~ 2 2 2 ~ ■ - ' “ 

- n k')(r - n k )A'(f) 

1 o 2 o 


where A (M / ) = e sin 2 ¥ + e cos 2 V . 

z 

Let the various parts of the z component of the field be defined 
as follows 

CO TT 

J^(r R sin e sin V ) 

(2.19) 


jn . l 2 f f . cos 2 ^ r 2 

ZS k o 2ir J j e sin 2 '? + e cos 2 '!' ° 

0 0 z 


-iTR cos 0 cos V . ... 

e J sin ¥ dYdr. 


. 2 r r (e 2 + e ,2 cos 2 '}')sin^ ¥ -jTR cos 0 cos ¥ 

E zfs = - jrlk o ( 2 rr ) J J A 2 ('F) 

0 0 


( 2 . 20 ) 


• J q (TR sin 0 sin V ) dTdr 

1 2 r .. r sin 2 '!' {f 2 k 2 [ e (e 2 -e ’ 2 )- (n 2 +n 2 ) (e 2 +£ * 2 cps 2 ^) ]+n 2 n 2 (e 2 +e » 2 cosV)k° } 

E zf = W'to? ] 2 : : 


0 0 


2 2 2 2 2 2 2 

(r - n O(r^ - n k ) A Z W 
i o 2 0 


( 2 . 21 ) 


t /nr, . „ . , -iTRcos 0 cos f . 

J Q (rR sm 0 sin V ) e sm V d'l'dr.. 


The total field is given by 


E = E + E £ + E 

z zs zfs zf 


( 2 . 22 ) 
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When the total field is decomposed in this way , E represents the most 

z s 

singular term in R, E - the next most singular term, and E' represents 

Z I S Z I 

a finite term. What remains now is the evaluation of each integral 
separately. 


2.3 Integrating E zg 


The convergence of (2.19), the equation which defines E , at first - 

zs 

appears to be questionable. One can note that for a three-dimensional 

2 

Fourier transform, V with respect to the observation coordinates R 

2 

and 0, can be replaced by -T . Therefore, (2.19) can be written as 

2 

= < 77 > V 2 I (R, 0) (2.23) 

zs k. 2 TT s 

o 


V l s (R, 0) = 



r 2 cos 2> y sin ¥ -jTR cos Y cos 0 


e sin 2 Y + e cos 2 Y 
z 


J Q (rR sin Y sin 0)drdY 
(2.24) 


where the order of integration and differentiation has been formally 
interchanged. 

Suppose we make the following change of variables to cylindrical 
coordinates. 


z = R cos 0 , 
y = T cos Y, 


p = R sin 0 
u = T sin Y . 


(2.25) 


Then we can write (2.24) 


oo Tf 


v 2 v r - e > 


cos Y/sin Y 

„ - e sin 2 Y + e cos“Y 
0 0 z 



e" JYZ J 
2,„ o 


dYdr (2.26) 


wher e Bess el’s equation Jias.Jaeen^ used -to- explicitly — p e r f o rm-d if f e r en tia t ion 
with respect to p outside the integral. 



13 


Completing the change of variable gives 


OO CO 


V 2 I s (p 


, U J. f f v 2 e' jYZ , , 

• z) p 3p P 3p J J ,2 2. J o 

..-n _ »(eu + e. Y > 


^ (up ) dydu'. 

n uCeu" + e y“) 

u =° Y~° * T (2.27) 

First, to perform the y integration, we will perform the integral 


OO 

■/ 


2 -jyz 




2 2 
eu + e y 
z 


dy. 


(2.28) 


At this point we can explicitly change the order of differentiation with respect 
to z and the integration,;, thus completing the formal interchange 
indicated earlier. 


J = 


-3 

3z' 


f e~^ Z 

J eu 2 + e y" 

-00 z 


dy. 


(2.29) 


At last we have an integration that can be performed with the aid of 
the residue theorem of complex variables. The integrand in Equation (2.29) 
is of exponential order; and therefore, for positive z we may close the 
integration path over the lower -half complex y plane. There is only 
one simple pole at 


Y = -j /e/e 2 u 


One obtains 




e zu 
z 


J = 


-3 tt e 


3z /e~ /e u 
z 


(2.30) 


We substitute into (2.27) and perform differentiation with respect to z 
-"Ve i a a °f 


V 2 I - 
s 


/s 


f_ i _3_ _3_ r 

3 p 3p P 3pJ € 


zu 


J (p u) du . 


(2.31) 


This final integration for E can be performed to obtain 

z s 
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-X = -1 $ JL 


' e z 3 p 9p "9 p J~2 


ez + e p 
z 


(2.32) 


where the following definite integral has been used [3] 

00 

..1 


e at J (bt) dt = 
o 


0 


7 2 ^, 2 . 

/ a + b 


(2.33) 


Finally 


= J 


ti Ve 19 9 


"zs 47rk^e_ p 9p P 9p J 2 


o z 


ez + e p 
z 


(2.34) 


or 


- 


zs 4irk e . 2 2 3/2 

o z (ez + e p ) 
z 


0 2 2 
2ez . - e p 
z 


2 , 2 
ez + e p 

Z 


(2.35) 


The dependence on distance can be seen more clearly by noting 


2 2 
5 +e p 
z 

; 

statics problems. Clearly 


that R' = /e /ez +e p is used as a plasma-scaled distance in some 
v z 


2 

E = 


zs 4irk e 
o z 


,.3 


, 2 2 
2ez - e p 
z 

i 2 , 2 

L ez + e p 
z 


(2.35a) 


2 . 4 Integrating E 


zf s 


We repeat (2.20), which defines E 


zfs ’ 


co 7T 


E _ = -jnk ( ) 

zfs J o 2 tt 



(e 2 + e’ 2 cos 2 y)sin\ -jTR cos 0 cos V 


0 0 


A 2 (y) 


J (rR sin 0 sin ¥ ) dTdV 
o 


AW' = e sin 2 T + e cos 2 T. 

z 
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This. term will be artificially divided into .two. :t-erms . an< * E zfs2 

for ease of integration. 

2 


E = -jnk ( ) £ 2 I- .. 

zfsl J o 2 tt fsl 


(2.36) 


where 


05 -Jf 


fsl 


'll 


0 0 


sin 1 -jTR cos 8 cos ¥ 

~2 e 

A^(¥) 


• J (TR sin 0 sin .¥ ) d¥dl;.. 
O 

E 2 £s2 ’ -J" k o < A > £ ' 2l £s2 


(2.37) 

(2.38) 


where 


00 "n" 


fs2 


’ If 


0 0 


cos 2 ¥ sin^ ¥ -jTR cos 0 cos ¥ T s 

o e J J (TR sin 0 sin ¥ ) 

A 2 (¥) 


(2.39) 

Again it will be useful to apply the same change of variables used for 


E . 
zs 


z = R cos 8 p = R sin 0 

y = T cos ¥ u = T sin ¥. 

Applying these changes, (2.37) becomes 

GO GO 

.3 


(2.40) 


u 


fsl 


u=0 , 


yL- oo (eu 4-ej ) 


2 2 ~ e " !YZ J C (P U ) dydu. 


(2.41) 


To do the y integration we need to evaluate 

CO 

-jyz 


-J-, 


2 , 2,2 
eu + e y ) 
z 


dy. 


(2.42) 


This integration can be done with the aid of the residue theorem. One 
can close the integration path in the lower -half complex y plane and 
evaluate the residue at the second-order pole, y = -j Jzft~ u, obtaining 


d¥dT. 
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J = 


-A 7 


e zu 
z 


2e e u 
z 


z + 


J e/e z u 


(2.43) 


Equation (2.41) becomes 


fsl 2e 


( U2 + 


i ] -As 


zu 




J (pu) du . 
o 


(2.44) 


or 


fsl 2e e 


z J_ 1 

- J e/t z 3z y e/e 


z J 


CO 

/ 

0 




e zu 

e J (pu) du 

o 


(2.45) 

by recognizing the integral as the derivative of a different function. 
Applying the integration formula 


00 

/ 


e aU J (pu) du = 1 

o 


/7 T ~ r 


(2.46) 


p + a 


obtained from Dwight [4], (2.45) may be written 


fsl 0 2 , 
2e / e 


" z 8l +1 


y-n i -i 

A p + e z 


-12 
z 
z 


(2.47) 


Performing the indicated differentiation, E ^ ^ ma Y be finally 


written as 


-jnk 


J zfsl 


, 2 . 2 . 1/2 
(e p + ez ) 
z 


P.+ 


2 

ez 


2 2 
e p + e z .J 
z J 


(2.48) 


or 


J zfsl 8tt R' 


1 + 


2 

ez 


2 2 
e p + ez . 
z 


(2.48a) 


The integration for E z f s 2 wil1 be done similarly. Using the same change 
of variables, 
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fs2 


r r 2 3 

'/ / <„*V. Y 2 ) 


u=0 y=— ' 

Again the y integration is done first 


o ^ e ~ |YZ J (pu) dydu; 
" / ™ “ \ o 

(u +Y )■ 


(2.49) 


r y2 e 3^ Z ( jY 

/ , 2 . 2,2. 2 . 2 ‘ 

J (eu + e y ) (u + y ) 


(2.50) 


This integration can be performed with the aid of the residue theorem 1 . 

Note in this case that there will be a pole at y = -u that did not 

occur for E - , . One obtains the following equation 
ztsl 


J = 


-IT 


e (e e)u~ 
z z 


f 

e e 
z 


uz 


e - e 

l 2 


+ e 


-/77 


£ UZ 

z 


e + e 

z 


2 /^>z " G) 


(2.51) 


The final integration of E , , uses the exact methods as E . > 

zfs2 zf si 

did. Differentiation with respect to z is moved outside the integral, 
the integration formula (2.46) is applied, and differentiation with 
respect to z is performed. After simplification one obtains 


jnk, t ,2 


zfs2 4tt (e - e) y (e - e) 


/ 2 2 2/e fe p 2 + ez* 

/z + p / z 


, _ 

e + e 

z 

e - e 2 2 . 

z e p + ez j 

7 


ez 


(2.52) 


or 


jnk ,2 


o e 


zfs2 4tt 


(e z’ E) \ (e z - e) 


/ 2 " 2' 2R ' 

/ Z + 0 


e + e 

z 

e - e 

z 


ez 


2 , 2 
e p + ez ■ 
z 


(2.52a) 
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2.5 Integrating; E 


zf 


Equation (2.21), which defines E r , can be written 

zl 


E rf - i" k o ‘ e > 


(2.53) 


IT 00 


P(R 


■•>-// 
0 0 


N(r ^ -i rR cos 6 cost. 

— S - 2 5 — 5 5 0 — ~ J (TR sin 0 sin ¥.) 

(r 2 - nf k 2 ) (r 2 - n 2 k 2 ) ° 

1 o 2 o 


• sin 3 T drd4“ 


(2.54) 


N(r, T) = r 2 [e(e 2 - e’ 2 ) - (n 2 + n 2 ) (e 2 + e’ 2 cos 2 Y)] 


L 2 2 , 2 . ,2 2.2 
+ n, n„ (e + e cos 4*) k 
l - L o 


(2.55) 


Note that the integrand of (2.54) is even in 4'. Therefore, the 
limits of integration can be changed 


P(R 


v, r 2 r N(r, ¥) sin 3 T:e jFRp J (TRq) 

, 6 ) = -5 0—0 j ° A ? dTd'i' (2.56) 

J 0 i A Z (4') (T 2 - n 2 k 2 )(r - n 2 k 2 ) 


where we have also substituted 


(2.57) 


p = cos 8 cos T 
q = sin 0 sin T 

The Bessel function can be changed to a more tractable form by using 
the integral definition 


J Q (TRq) 


TT / 2 

= ’ l 


a — j TRq cos a + ^jTRq cos a 


da. (2.58) 


Making this substitution and changing the order of integrations, we 
can organize P(R, 0) as 
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where 


tt/2 u/2 


P«,e)=i J J( i 1 + i 2 ) 


0 0 


sin 3 y 
2 ' 

A Z (T) 


dady 


(2.59) 


X 1 = 


CO 

I 


N(r, ) e 


-j TR(p + q cos a) 


, r 2 2 . 2 2 2 . 2 . 

-co < r " n i k 0 > < T ~ n 2 k o ) 


dr 


(2.60) 


V 


N(r t f).e- jTB(p ~ q COS a) dr< 

/p 2 2.2 Wp 2 2,2. 

(F - n r k ) (l - n_k ) 

1 o L o 


(2.61) 


The next step is obviously integration by the residue method . 

However, one must take care to consider the signs of (p + q cos a) 

and (p - q cos a) for various ranges of ¥ and 6. Also one must 

2 2 

remember that the square roots of n^ and n^ must be taken so that n^ 
and n 2 have negative imaginary parts to satisfy the radiation condition. 

For the range of the angles; permitted by (2.59) for I^,it can 
easily be seen p + q cos a > 0. Therefore, the integration path can 
be closed in the lower-half T plane to obtain 

-jn k R(p + q cos a) -jn k R(p + 

N (n k ,T) e N (n k ,T) e 

1 o , l o 

2 2 2 2 
2n 1 (n 1 - n 2 ) ^n^^ - n 2 ) 

(2.62) 


i = _ 2 ii 
1 k 3 


For I 2 , the sign of (p - q cos a) depends on the exact range of 
0, f, and a. Recalling (2.57) 


p = cos 0 cos ¥ 
q = sin 0 sin ¥, 


(2.57) 


cos 


one can see 
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p - q cos a > 0' whenever either 0 < 4* < tt/2 - 9 and"0 < a < tt / 2 or whenever 


tt / 2 - 9 < ¥ < tt / 2 i and cos ^ ^ <a< n/2. 

- q - - 


(2.63) 


In this case the integration path can also be closed in the lower-half F 
plane to obtain 


I = 

2 k 3 
o 


N(n,k , « - jn l k „ R(p - 1 cosa) N(n,k , T) -3 n 2 k „ E(p “ 1 cos “> 
1 o 1 ° 


1 o’ 


2 "l<-\ - p 2> 


-2n 2 (n 2 - n 2 ) 


(2.64) 


However, . ' < •- 

p - q cos a < 0 for - 0 < V < ^ and 0 < a < cos . 

(2.65) 

Then, the integration path must be taken in the upper half of the complex - 
T plane to obtain 


I = 2rri 
2 k 3 


O L_ 


jn..k R(p - q cos a) jn k R(p - q cos a)-, 

N(- ni k , V) e 1 ° N(-n„k , V) e 

1 o . l o 


-2n 1 (n 2 - n 2 ) 


2n 2 (n l ” n 2> 


( 2 . 66 ) 

Gathering the results of (2.62), (2.64), and (2.66) and carefully 
noting over which range each is valid, we can rewrite (2.59) . Note that 
N is an even function of its first argument, and that the second term of 
1^ and I2 is the first term with all n^'s and ^'s interchanged. 


P(R, 9) = R nl (R, 9) + P n2 ( R » e > 


(2.67) 


P „(R, 0) = P , (R, 9) with n, and n„ interchanged 
n2 ni 1 1 


tt/2 tt/2 


f nl «. e> 


i>-* f [ 

k J. J 


N(n 1 k Q , 4') e 


-jn^k Q R(p + q cos a) 


_ n l (n l ’ “2? 


sin 3 f 

A 2 W 


dadV 


21 


«<n, k , « e' 30 ^ ' q “* B) , 3 

v n „> ' sm _xjr 


. j . r r mu iV *- 

k 3 X 1 n 


o 0 "0 


2 2 

l (n l - n 2> 


A 2 (¥) 


dadY 


tt/2 tt/2 N(n.k , Y) e 

If 

Ar, „ J -1 


-j n lk o R(p - q cos a) 


o tt/2-8 cos p/q 


2 2 
n-j 0^ - n 2 ) 


• 3„ 

s in V 

a 2 (y) 


dadY 


it/ 2 


, / / 

o tt/2-0 0 


-1 , +jn,k R(p - q cos a) 

c °, s P'" N(n l *) e 1 ° 


, 2 2 , 

n 1 U 1 - n 2 > 


sin~^Y 
A 2 (4?) 


dadY, 


( 2 . 68 ) 


Note that the first integral* s range can be broken into the ranges 
of the last three integrals. One can then distribute the integrand 
of the first integral to the last three. Recognizing the exponential 
form of the cosine function one finds 


P nl (R ’ 6) = 


tt/2-0 tt/2 


0 ^0 


" jn l k o Rp 

H(n^, n 2> 40 e cos(n^k_Rq cos a) dadY 


tt/2 tt/2 


tt/2-0 cos ‘'‘p/q 


H(n^, n 2 , Y) e 


-j n l k o R P 


cos(n v k Rq cos a) dadY 
1' o ^ 


tt/2 cos '‘p/q 

+ J f n 2 , Y) e 

* — - y . 


-jn k Rq..c°s a 


tt/ 2-0 


0 


cos(n,kRp) dadY 
I o 

(2.69) 


where we have let 
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N(n T k , 40 sin ¥ 

H(n , n , 40 = ^ — — 

x z k J 2 2 2 

o A^(40 - np 


(2.70) 


A little algebra simplifies (2.70) to 


H(n^, n 2> ¥) = k 


r , 2 , 2 . 2 . 2 , ,2 2 . . . 

-2‘ n ]0 £ ( £ - £ 7 — n^(e + e cos 4')]sx 


sin^y 


2 2 2 
A 140 (n^ - np 


(2.70a) 


Now if we formally add and subtract the term 
tt/2 


cos lp/q 


j f H(n 1 , n 

tt/2-0 0 


- jn l k o Rp 


2 » 4*) e cos(n^k Q Rq cos a) dad4' 


from p n ^( R > 0) , it completes the range of one integral so that we can use the 


identity [5] 


it/ 2 


V"i 


k Rq) - - f 

O. IT / 

7n 


cos(n^k Q Rq cos a) da. (2.71) 


Using these two manipulations jP n ^(R, 0) further simplifies to 

tt/2 


P ni (R . 0 ) = \ J H(n lf n 2 , 4-) e 
0 


"j n l k o R P 


J (n 1 k Rq) d¥ 
o 1 o 


tt/2 cos ^p/q 


J J H(n 1 ,n 2 ,'R) [ 


-jn^k^Rq cos a 


cos(n^k Q Rp) - e 


-jn..k Rp 
l.o 


tt/2-0 0 


cos(n^k Q Rq cos a)] dad4'. 


(2.72) 


Simply expanding exponentials 

tt/2 

P nl (R ' S, ’T { H <V n 2’ n e 


-Jn^Rp 


J (n..k Rq) d4' 
o 1 o ^ 


tt/2 cos 1 p/q. 


+ j 


H(n^, n 2> 4') sin n 1 k Q R(p - q cos a) dad^.. (2.73) 


tt/2-0 


0 
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At this point it appears that no more exact simplifications are 
possible. It will be necessary then to make judicious approximations to 
put (R, 0) into a suitable form for calculation. The approximations 
will be easier to justify in cylindrical coordinates where 


z = R cos 0 
p = R sin 0. 


(2.74) 


Recall that z is to be interpreted as the distance along the antenna - . 
from the dipole source and p will be a, the radius of the antenna, which 
is assumed to be small compared to the length. 

First, note that the sine term in the second integral of (2.73) is 


sin n 1 k [z cos V - a sin *P cos a] = sin (n.k z cos V ) 

1 o 1 o 

* cos (n, k a sin *P cos a) - cos (n.k z cos *P) sin (n,k a sin I cos a).. 
1 o 1° 1 o 

(2.75) 

For small a we make a first-order approximation of (2.75) by 

sin (n k z cos '?) - n.k. a sin. 'P cos a cos (n.k z cos . *P) , (2.76) 

to l o i o 


Now we can perform the indicated, a integration in (2.73) on this simplified 
term (2.76) to get 


z cos 'F 

sin (n r k z cos *P) cos — ■ — : — nr - n.k a sin 'P cos (n,k z cos *P) sin a 

1 o a sin *P 1 o 1 o 


z cos ¥ 


where / 

sin a = /l - 

/ a sin 'P 

We shall also make an approximation for the Bessel function in (2.73) 


(2.77) 


J (n.k Rq) = J (n,k a sin f) =r 1 
o lo o lo 

for small a. Finally, we can write 


(2.78) 
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P nl (R ’ 0) 


-jn,k z cos ¥ 
1 o 


V 2 

/ H(n^, n^, W) sin(n^k o z) cos ¥ cos ^ 

tt/2-0 *• 


-1 z cos y 
a sin ¥ 


n,k a sin y cos(n..k z cos 4*) /l - z c ? s ■-] d'i' (2.79) 

1 o . 1 o / a sin ¥ 


where 


. -la 

9 = tan — . 

z 


2.6 Final Field Expressions for Infinitesimal Source 


At last all the fields from an infinitesimal source have been found. 
They are repeated here for references purposes. 


E = E + E ' ' + E . 0 + E . 

z zs zfsl zfs2 zf 


(2.80) 


4rk e D , 3 

o z R 


, 2 2 
2ez - e p 
z 

2 , 2 
ez + e p 
■ z 


(2.81) 


J zfsl 8r 


£ . -L • 


2 , 2 
ez + e p 
z 


(2.82) 


5 ^ ,2 r i 

"zfs 2 4?r e-e 1 e - e 
z z 


_ | r / 2 2 

R = /e J ez + e^p 


2 2 
s + p 


(2.83) 

(2.83a) 
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where 


H(n 1 , n 


E zf ■ « k o < >' l p nl (z - "> + p „2<“- p)] (2 - 84 > 


o ' 2ir ' ‘ nl 
it/2 


P nl <2 > p) ’ 


-* / 

rv 


H(n 1 , n 2> 40 e 


-jn^k z cos 4* 


d4* 


+ J 


1t/2 

/ 

ir/2-8 


H(n 1 , n 2 , 40 


. / , ,„S -1 z COS 4 1 

sm(n,k z cos 40 cos 7— — 

1 o a sin 4* 


n,k a sin V cos(n_k z cos 40-\/l 
1 o 1 o 



z cos 4* 


a sin 4' 


d4 f 


(2o 85) 


, V) 



6 = tan ^ — 
z 

r t 2 » 2 \ 2/ 2 . .2 2 . . 3... 

[e(e - e ) - n^(e + e cos 4')] sin f 

2 2 2 
A Z (4')(n^ - np 



2 B(4Q 
n 2 



(4Q - 4A(4QC(4Q 
2A(40 


(2.85a) 


; Im(n 1 ) < 0 , Im(n 2 ) < 0 

( 2 . 86 ) 


2 2 

A(40 = e sin 4' + cos 4' 

B(4') = ee (1 + cos 2 4') + (e 2 - e 1 2 ) sin 2 4* (2.87) 

25 

C(4*) = (e 2 - o' 2 ) e (2.88) 

z 

P n2 = P n ^ with all n^'s and n^' s interchanged. 
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CHAPTER III. APPLICATION OF METHOD OF MOMENTS 

In the last section we derived an expression for the electric field 
produced by an infinitesimal source. It can easily be seen that the 
current distribution on the antenna, i(z), must satisfy the integral 
equation 

i 

i(z’) K(z - z') dz' = e(z) (3.1) 

o 

-E applied at antenna feed, 

where e(z) is 

Lq everywhere else 

apd where i is the length of the antenna, and K(z - z') = E (z - z') is 

z 

the Green's function evaluated on the surface of the antenna. Thus 
(3.1) basically specifies that the tangential E field on the surface 
of the antenna must be zero. 

' i . ' 

The general theory of the method of moments will be presented, following 
the formulation by Harrington [6]. Then it will be shown how this method can be 
applied to this problem. First define the inner product of two functions on 
some interval I as 

Now consider solving the inhomogeneous equation 

' i . 

L(f) = g (3.3) 

Where f and g are functions defined on the common interval I, L is 
a linear operator, g is known, and f is unknoWn. 

Suppose f is approximated _as_ a Siam of n known_ f unc ti on s _we i gh t e d_b y — 

n unknown coefficients a . 

» n 


< f ( z) , g(z) > 


f(z) g(z) dz. 
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£ = Z cj V 
j = l 


(3.4) 


The fj functions are called expansion functions. Now by substituting 
(3.4) into (3.3) and using the linearity of L, the equation becomes 

(3.5) 


n 


l « L(f ) = g. 
j=l 2 2 

Now define a set of n weighting functions or testing functions, 
w^, i = 1, n. By taking the inner product of both sides of (3.5) with 
respect to each testing function, the equation becomes 
n 

I a < w , Lf > =? < w ,g>, i «• 1, 2, • • • n. (3.6) 

j=l * 1 2 

This can be written in matrix form as 


UijHpij] = tg ± ] 


(3.7) 


where 


[ V ■ 


< w x , Lf x » < w x , Lf 2 > 


< w 2 , Lf x > < w 2 , Lf 2 > 


t • • 


[0j ] 


“l 


8 i 


< w 1 , g > 

a 2 

o 


g 2 

• 


< w 2 , g > 

• 

[ gi ] = 

• 

B8 

c 

a 

n 




« 




J V 8 L 


(3,8) 


(3.9) 


Equation (3.7) can be solved for [ a^ ] simply by inversion, provided 
[£ '] is. nonsingular 

ij 
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ta j ] = (3.10) 

Thus, f has been solved as a weighted sum of the expansion functions. 

The method of moments, as it has thus been described, is very general. 
The particular choice of expansion and testing functions used for this 
problem will simplify the form of the solution considerably. Also, 
this choice affects the accuracy attainable for a given number of 
functions. Now we will show how to; apply the method to this antenna 

problem by first identifying the functions and operators in (3.3). The 

i 

unknown function f is i(z) defined on the interval [0, £.] representing 
the antenna. Likewise, the known function g is the tangential electric 
field e(z) . The linear operator L is the integral operator 

a 

L(h) = J h(z’) E z (z - z') dz' (3.11) 

o 


or described in words, convolution with the Green's function along the 
antenna. 

Next expansion functions must be devised for i(z). We will follow 

what is called the method of subsections, that is, expansion functions 

are chosen that are defined on only a subsection of the antenna interval. 

First, the antenna is divided into odd n number of segments, where each 

2ti 

segment has length — . Then, the expansion functions are defined as the 


pulse functions 


'1 if (i - 1) — < z < 1 ( — ) 
n - n 


P^z) 3 ' 


. (3.12) 


Q otherwise 


Now we can express i(z)_as 


i(z) 


n 

I 

i=l 




(3.13) 
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Figure 2 shows how the use of pulse expansion functions replaces the actual 
current with a step approximation. 

Now a set of testing functions must be chosen. The simplest choice 
utilizes the point-matching method. For this method, we choose w i to 
be 6(z - m^) , where 6 (z) is the Dirac delta function and is the 

midpoint of the i 1 "* 1 antenna segment. Thus, by the sampling property 
of the 6 function, the right-hand side of (3.6) is 


< g > = gCmp. 


(3.14) 


Therefore, for this problem [g^] = [e^], where [e^] is the vector of 

tangential electric fields evaluated at the midpoint of each segment. 

Similarly, we can see that element of (3.8) is the electric field 

tih 

at the midpoint of the i segment caused by the uniform current over the 
.th 

j segment. 

Finally, the original integral equation (3.1) has been transformed 
into the following matrix equation 


Z I 


(3.15) 


where Z = [z^]j z ij t * ie electric field produced by the uniform current 

til til 

over the j segment evaluated at the midpoint of the i segment. 

0 

0 

0 


I = 


i, is the current 
k 

i th 

over k segment 


E applied 
0 
0 
0 


tangential 

electric field 

at midpoint of 
, th 

k segment 


The solution is obviously 


I = Z -1 E = YE. 


(3.16) 



interval 


2 


PULSE FUNCTION Pi (z) 



2 


PULSE FUNCTION EXPANSION FOR i(z) 

Figure 2. Use of Pulse Expansion Functions 
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Suppose the feed is at the k segment. Then the impedance of the 


antenna is simply - — where v, is the voltage over the gap. In this 

\ k 

report the feed is always at the center of the antenna. 


Note some very important properties of the matrix Z. First, there 


are only n unique elements. Segments the same distance apart have the 


same effect on each other. Thus z^ = z^ + ^ j+k' Also, by symmetry 

z., = z... Therefore we can write Z as 
ij Ji 



Z = 


. ( 3 . 17 ) 


z 

n 



The matrix Z is a complex symmetric matrix. Unfortunately, most 
of the useful properties of real symmetric matrices or complex 
Hermitian matrices do not apply to complex symmetric matrices. The 
element z^ is called the self-element since it is the field at the 
midpoint of a segment carrying a uniform current. Since the fields 
decrease strongly with distance, it will be shown that Z is strongly 
diagonal. 

The Z matrix has been discussed somewhat generally. Next, we 
will discuss the method for calculating the elements from the Green's 
function derived in the last section. There, are three different sets of 
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expressions based on three different models for the distribution of 
current over one segment. The basic model used assumes that the 
current over one segment is uniformly distributed over the axis line; 
this current distribution will be called the current filament model. 

Since the field is evaluated at the radius, it remains finite. 

Recalling that the length of one segment is 2A£ , we can write the 
field caused by the current filament as the superposition of fields caused 
by infinitesimal sources distributed along the filament. Thus, 


E cj i _(z, a) 
z filament 



z' , a) dz' . 


(3.18) 


The factor of 1/2A£ is included to make the effective current over the 
segment equal for the two models . This integration was done on each 
of the four component fields of E^. For E zg , E^^, and ® z f s 2» one 
must perform integrals of the form 


A£ 


■L 


(z - z')" 
t(z - z') 2 + a 2 ] m 


dz 1 


(3.19) 


where n is 0 or 2 and m is 1/2, 3/2, or 5/2. These integrals can be 
found in common integral tables. Therefore, the integrated forms of 
E , E f ^ , and E f „ will be simply presented later. 

ZS ZIS _L ZXS / 

Integrating E ., the component of the Green's function E , which 
zx z 

cannot be written in closed form, requires more explanation. Recall 
that E^ is given by (2.84) 

E zl = Jnk o ( 2? > If nl Cz > p> + p n2 (z ' p>) - 

(2.84) 

The first term can be written as 
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P , = I , + I _ 
nl zfa zfb 


(3. 20) 


where 


zfa 

ir/2 


tt/2 

-* L 


-jn.k z cos ¥ 

HO^, n 2 , ¥) e 1 ° dr (3.21) 


zfb 


= j J H( n;L , n 2 , V) sinCn^z cos ») cos" 1 ) f ff^~f [ 


where 


tt/2-0 

- n.k a sin ¥ cos(n..k z cos 

1 O 1 0 


9 = tan 1 — . 

z 




df 

(3.22) 


Remember n^, n 2 and H are functions of f » 

We will make the approximation that I^.- can be neglected „ To 
justify this, first note that the integration width for the self- 
element where z = 0 is 9 = |, For the first segment away from the j . 
source element, z = 2A£, and since only thin wires are being treated, 
— is small. Therefore for all elements, except the self-element s 6 - 
We can now make small angle approximations in (3.22) and perform the 
integration to get 

2 


I zfb ~ "I H(n l* n 2* 2^ 8 n l k o z » 


z > 2A£. (3.23) 


' SL 

Because of the — term, I _■ should be small compared to I f . For 

; Z ZXD Z X cl 

z =0, Equation (3.22) simplifies to 

ir/2 


X zfb " j 


H(n^, n 2> 1/) n^ k Q a sin Y dr (3.24) 


■ , i ; 

Direct numerical integration of (3.21), (3.22), and (3.24) indicates 
that I ^ may he neglected for all terms including the self-element. 


Ctf.| N 
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Moreover, one should note that the singular field terms are orders of 
magnitude larger than E ^ for all but the farthest segments; therefore, 
E^ need not have an extremely small relative error. 

Now, when I ^ is neglected, E^ can be integrated over one 


segment, according to (3.18). Note that z occurs only in the exponential 
-jn^k z cos 'F 

e . Therefore, one obtains the integrated form of E^ 

by integrating the exponential. One can easily verify that 

fsinCn^k^cos 'F AJl) 


A £ 

-i- f 

2A£ J 
-ML 


-1n..k (z - z') cos ¥ -In.k z cos *F 

J 1 o , , J 1 o 

e dz’ = e 


. n n k cos 'F AJl 
L 1 o 


(3.25) 

One does the exact same integration for the exponential containing 
n^. The integrated form of E^ is then only E^ with the exponentials 
replaced by thh right-hahd side of (3.25). 

We will summarize the fields produced by a current filament mode. 
To avoid proliferation of variable names, the integrated form has 
the same name as the corresponding term of the Green's function. We 
have substituted the antenna radius a for p. 


z filament model zs zfsl zfs2 zf 


(3.26) 


JiL 


z ML 


z - ML : 


J zs 8Trk o eAS, [ [ (z + AJl) 2 + a ,2 ] 3/2 [ (z - AJl) 2 + a' 2 ] 3/2 

(3.27) 

,2 e z 2 
l = 


where a = — a 

e 


7 = ijV 

J zfsl 8tt (2A£) 


2 Jin 


+ ML + /(z + Ail) 2 + a' 


z - ML 


♦ / 


(z - Ail) 2 + a' 2 


z + AJl 


z - ML 


J (z + 


AJl) 2 + a' 2 


/ (z - AJl) 2 + a' 2 


(3.28) 
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J k o n £* 2 1 


/ 2 
+ /(z + A£) + 


J zfs2 4 tt (e - . e) 2A£ ye - e 


JL ^ z + ,.A£ + /(z 4 A£) + a 


Z ' Z z - A£ + J (z - A£) 2 + a 2 


1 „ ' z + A£ + /(z +A£) 2 + a’ 2 

£n 

+ /(z - A£) a’ 


e -- e A „ 

z z - A£ 


1 

2e 


z + A£ 


z - A£ 


/(z + A£) 2 + a ' 2 /(z - 


(3.29) 


(z - A£) 2 + a' 2 


E zf = ^ k o ( 27 )2 (P - + P - } 


L nl n2‘ 


(3.30) 


nl 


ir/2 . . 

r “jn.k z cos ¥ 

~ j H( ni , n 2 , >F) e ° 

J o 


sinOyk cos ¥ A£) 
1 o _ 

n^k cos ••¥ A£ 

1 o 


d'F 


(3.31) 


H ( n l* n 2’ ^ “ k ' n l 

o 


r/2 ,2. 2 2 ,2 2 .... . . 3 

[e(e - e 1 ) -r n^(e + e cos ¥)] sin 'F 


2 2 2 
A (¥) (n x - n 2 ) 


(3.32) 


i±/e 2. 


n 2 , 2„ BCf)*/-B ; CQ -4AW C(f) . < 0> (3.33) 

-■ 1 2A('F) X 


Im(n 2 ) < 0 


A(¥) = e sin 2 ¥ + e cos 2 ¥ 

z 


(3.34) 


B(¥) = ee (1 + cos 2 ¥) + (e 2 - e’ 2 ) sin 2 ¥ 
z 


C(¥) = (e 2 - e' 2 ) e 


P _ = P., with all n 's and n's interchanged. 
n2 nl 1 2 


(3.35) 
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Equations (3.32) to (3.35) are the same as Equations (2; 80) to 
(2.89) but are reproduced for reference. 

It has been shown [7] that one obtains better answers if /2a 
is used in the calculations as an effective radius in place of a. 
This correction compensates for the assumption that the current is 
on the axis instead of distributed on the surface. 


One can note that many simplifications in the preceding : formulas 
are possible and in many cases necessary for limiting cases of the 
general anisotropic plasma. For free space, the E ^ 2 term has zeros 
in both the numerator and the denominator. One can clearly see that 


for plasma parameters X = 0, corresponding physically to free space 

2 

with an applied magnetic field, and that for Y =0, corresponding 
physically to an isotropic plasma, there is no off-diagonal dielectric 
element e' and, therefore, = 0. In the limiting case of 

free space, E^^ must be set equal to zero. One must also analytically 
treat removable singularities in the expressions for E^. For the 
isotropic case, the two indices of refraction n^ and are equal. 

The expression for H can be rewritten to avoid an indeterminate form. 

Next one can compare the fields produced by the current filament- 
model to the fields produced by 6-source model. The 6-source model 
is physically unrealistic; however, it should give the same values, as 
the filament model for observation points far from the segment 
with the current. Far from the source, numerical calculations show 
that for all field components one may use the 6-source model to avoid 
subtracting nearly equal terms present in the filament model expression. 

^In_general,_no_more_than_21_segments-were used -and,— therefore,— the 

6-source model was not needed for computation. One can note that the 
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E , term is nearly independent of distance. Therefore, one would 
zf 

not expect the fields due to the two models to be significantly 

different for any distance. Calculations verify this. 

A third current distribution model was considered. One might 

object to the filament model as being unrealistic, especially for 

elements very close to the current segment. Therefore, we have 

examined a model for which the current is uniform on the surface 

of the segment. We call this model the cylindrical current model. 

This model can be obtained from the filament model by summing up the 

contribution of distributed current filaments placed along the surface 

of the current segment. When one integrates around the segment, one 

needs to replace the radius a with the radial distance 2a sin 

(See Figure 3.) Recall that E has the form E (z, a). Therefore, 

z z 

we write 

2tr 

E , . . . , = ir~ I E (z, 2a sin ^-) dtp . 

z cylindrical 2ir I z filament 2 

0 (3.36) 


Since the cylindrical mpdel is needed only in the very near region, and 

since in this region the E term strongly dominates the other terms, 

z s 

this cylindrical integration is applied only to the E term. Then 

z s 

substituting E^ f ^ 1 ampr)1 _ into Equation (3.36) one obtains 


2t r 


'zs cylindrical 2tt 
z + AZ 




8frk eAZ 
o 


z’ - AZ 


L [ (z + A£) 2 + 4a ,2 sin 2 


3/2 


[(z - AZ) 2 +• 4a ,2 sin 2 |] 


3/2 


dp 


(3.37) 
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PROJECTION INTO x-y 


Fi gure 3 . C ylindrica l Model as Pis tribute d _Cur rent _ Filaments 
of Segment 


PLANE 


on^Surface 
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One can see that (3.37) essentially is an integral of the form 


I 


J_ 2 ? <a« 

2 * J 0 [4a- 2 sin 2 f +b 2 ] 3/2 


Making the change of variable 
and simplifying yields 


I = 


-J—k 3 f n 

4713,3 J 0 


6 = $ + IT 


de 


(i - 


72 . 2 .3/2 

k sm 0) 


(3.38) 


(3.39) 

(3.40) 


where 

k = — 9 2a 'o -,/ T and a ,2 = — a 2 . (3.41) 

0> + 4a ' )' 

Note how similar (3.40) is to the definition of an elliptic 
integral. Applying a clever transform from Gradshteyn and Ryzhik 
[8] allows us to write . 


I = — ~ k 9 E(k) (3.42) 

4na ,J 1 - k 

where E(k) is an elliptic integral of the second kind. Equation (3.42) 
can be applied to (3.37) to obtain 


zs cylindrical 
k 3 

(z-AS) - 

3 2 

4Tra' 1 - k 


_JH: 

8nk eAZ 
o 


z + A& _ 
4ira' 3 1 


E(k_) 



E(k + ) 


(3.43) 


. = 2a^ 2a 1 

K, 2 21/2’- 2 21/2* 

[ (z + AZ) + -bfi'T. [(z - MLV + 4a’ ] i/Z 

Now one can compare the fields due to a current filament to the fields 
due to the equivalent current flowing on the cylindrical surface. 
Although one might expect much different values for the most singular 
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terms* the results were extremely close. Except near z = Ail, the fields 

generally agree to better than 1 per cent. See Table 1 for a comparison 

of the two models for various distances and plasma parameters. 

One can note that the field for the cylindrical model is singular 

for z = A£, and that it changes sign crossing this point. This can 

easily be understood from a quasi-static model. From this viewpoint, 

the current causes rings of charge to accumulate on the edges of the 

cylinder, and the electric field is the electrostatic field due to these 

two rings. Clearly for z = AJl, one is evaluating the field at the 

ring and it is singular there. The filament model remains finite at 

z = A£ because for this model the charge is on the axis, an effective 

distance of a' from the point of observation. 

After comparing the three current distribution models, we 

decided that the filament model was entirely adequate for E ^ 

E , and E ... For E we use the cylindrical model although it 
zr S 4 zr zs 

is only slightly different from the filament model. Therefore, 
each element z^ of Z, as shown in (3.17), is computed as 

z . = E [ (i - 1) 2A£ , a] (3.44) 

1 z ■ . 

where 

E = E (z, a) ; 
z z 

the cylindrical model of E is used for E and the filament model 

z zs 

is used for E^, E^, and E if . 

There are still a few small matters to discuss before calculating 
I, the current vector. First, one may consider applying one volt to 
the-f eed obtaining -the- resulting -electric -field : 
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TABLE 1. 


COMPARISON OF ELECTRIC FIELDS DUE TO FILAMENT AND 
CYLINDRICAL CURRENT DISTRIBUTION MODELS 

H = 1, a = 1/144., A£ = yj k Q = .1, z = m2A£ 


2 

X Z 

segments 

m 

E cylindrical (z, a) 

z s 

E 

zs 

filament (z, /2 a ) 



REAL 

IMAG 

REAL 

IMAG 

0 0.0 

0 

0.0 

. 3395xl0 7 

0.0 

. 3388xl0 7 


.43 

0.0 

.2737xl0 8 

0.0 

.1418xl0 8 


1 

0.0 

-.1496xl0 7 

0.0 

-.1493xl0 7 


2 

0.0 

-.1283xl0 6 

0.0 

-.1283xl0 6 


5 

0.0 

-7437. 

0.0 

-7437. 

/ 

10 

0.0 

-916.8 

0.0 

-916.8 

•1 .1 .Of 

0 

-.2862x10"* 

. 3814xl0 7 

-.2835xl0 5 

. 3806xl0 7 


.43 

-.2067xl0 6 

. 3059xl0 8 

-.8602xl0 5 

.1574xl0 8 


1 

.1260xl0 5 

-. 1681xl0 7 

.1246xl0 5 

-.1677xl0 7 


2 

1094. . 

-,1442xl0 6 

1093. 

-. 1442xl0 6 


5 

63.54 

-8362, 

63.53 

-8362. 


10 


7,834 


-1031 


7.834 


-1031 
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TABLE 1„ (CONTINUED) 


2 

X Z 

segments 

m 

E cylindrical (z,a) 

zs 

E filament (z, /l a) 
zs 



REAL 

IMAG 

REAL 

IMAG 

.9 .9 .05 

0 

-.4836x10^ 

-,4101xl0 6 

-.4838xl0 6 

-.4099xl0 6 


.43 

1258xl0 8 

-. 9702xl0 7 

- . 1283x10 8 

~.9433xl0 7 


1 

. 2150xl0 6 

. 1822xl0 6 

.2150xl0 6 

.1821xl0 6 


2 

.1718xl0 5 

. 1459xl0 5 

,1718xl0 5 

.1459xl0 5 


5 

} 

985.6 

837.4 

985.6 

837.4 


10 

121.4 

103.1 

121.4 

103.1 

.2 .25 .0! 

0 

i 

-. 2767xl0 6 

-,2084xl0 7 

-.2754xl0 6 

-.2082xl0 ? 


.43 

r .2534xl0 7 

- . 2088xl0 8 

-1.440xl0 7 

-,1382xl0 8 


1 

. 1223xl0 6 

£ 

.9215x10 

.1216xl0 6 

,9207xl0 6 


2 

.1028xl0 5 

.7692xl0 5 

.1027xl0 5 

. 7692xl0 5 


5 

594.0 " 

4441. 

594.0 

4441. 


10 


73.20 


574,2 


73.20 


547.2 
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TABLE 1. (CONTINUED) 


segments 

m 

E cylindrical (z, a) 

z s 

E filament (z, /2 a) 
zs 

• 


1 

IMAG 

REAL 

IMAG 

0 

-.1291xl0 6 

. 1250xl0 7 

-.1312xl0 6 

. 1249xl0 7 ' 

.43 

-,.1857xl0 8 

. 2765xl0 7 

-. 2355xl0 9 

. 8449xl0 8 

1 

i. 

.5776xl0 5 

5571xl0 6 

.5877xl0 5 

-. 5568xl0 6 

2 

l 4361 - 

-4350xl0 5 

4372. 

-. 4349xl0 5 

5 

1 

; 248.4 

| 

-2488. 

248.5 

-2488. 

10 

|30.55 

i 

\ 

-306.2 

30.56 

■ 1 

-306.2 

1 

0 

i 

^•.1180xl0 6 

.2381xl0 6 

-.1180xl0 6 

t 

.2 380x10 7 

.43 

-.1499xl0 7 

.3021xl0 8 

-.1208xl0 7 

.243xl0 8 

1 

.5231xl0 5 

-. 1055xl0 7 

.5231xl0 5 

-. 1055xl0 7 

2 

,4288. 

-.8647xl0 5 

4288. 

-. 8647xl0 5 

5 

426.9 

-4979. 

246.9 

-4979. 

10 

30.41 

-613.4 

30.41 

-613.4 



















since the gap distance d is 2A£. Second, note that E f-Q ament gi ves 
the field due to a distributed source with the same effective current 


as the 6-source model o This 6 source was one ampere at one point; 
when it is distributed, the current at any point is only ampere. 
Therefore, one should multiply all elements of Z by 2A£ to get true 
fields. Equivalently, this last factor of 2A£ is factored from 
the Z matrix into the E vector. Therefore, considering both comments 


about the 2 A & factor, the vector E is computed as 



i = center element 
i ^ center element 


(3.46) 


One must also make slight corrections between effective and 
actual dimensions of the antenna. The use of pulse expansion functions 
assumes a finite current at the end of the antenna when the true current 
must be zero. To compensate for this, dummy segments have been 
added to each end. These segments are not used in calculations, 
but only make the length of the antenna for which the nonzero currents 
can exist shorter than. the actual physical length 2H. Therefore, -for 
segments one uses 


2h 

instead of 2A& = — . 

n 


The exact number of segments to add was partially determined 
empirically. Using one extra segment on each end results in a 
good, but not perfect, linear current _distribution for. a_very_ short 
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antenna in free space. The other dimensional correction which 
has already been noted, is that one uses an effective radius of 
1 Fi a instead of the actual radius a in all calculations using the 
filament model. 

In most cases,. 21 active antenna segments were used. Trials 
with 33 and 45 segments show only slight changes. It is not at 
all clear what is an optimum number of segments to use because one 
cannot always assume that using more terms gives a better approximation. 
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CHAPTER IV. NUMERICAL DIFFICULTIES 

The field expressions for the filament current model appear 
deceptively simple. Even though these expressions have required a . . 
lengthy derivation, most of the effort expended on this problem has been . 
in performing numerical computations. For the terms in closed form, one 
must guard against loss in accuracy due to subtracting nearly equal terms,. 

and for the terms indicated as integrals, one must perform difficult 

numerical integrations. The purpose of this section is to describe the 
numerical analysis required in writing an efficient computer program. 
Particular attention is paid to the numerical integrations needed for 
E^. The general question of stability of the method is discussed in 
Chapter V. 

The computer program that is used can be found in the Appendix. 

For computational purposes, it was convenient to formulate the field 

expressions with dielectric constants normalized to e. Thus, in the - 

program e' is actually e'/e, e is e /e, and the field expressions have 

z z 

been correspondingly modified. One result of this renormalization is 

that one uses n^/e in place of n^. In this case, one must be careful 

2 

to compute n^/e as the square root of n^e with the negative imaginary. . 

2 . 

part rather than as the product of the square roots of n^ and e taken 
individually. 

For the singular terms in closed form, E zg , E z f s ]_> and E^^, one 
must be careful to avoid subtracting nearly equal quantities. This problem 
increases with increasing distance z: As has been already indicated, if 

fields are needed more .than 30 segments away, the 6-source model gives 
entirely adequate values. A series expansion of the singular terms has 
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been made, but the range of validity for the filament and the 6-.source 
models overlap and, therefore, the series expansion is not needed. 

The largest .source of numerical difficulty is performing. the 

numerical integration for the finite term E^. The shape of the 

graph of the integrand varies with the three plasma parameters,, the .free- 

space wavenumber k Q , and the distance z. Figure 4 shows real and imaginary 

parts of the integrand for various parameters scaled so that the difference 

between maximum and minimum function values is unity. For some plasma 

parameters the integrand is well-behaved, but for many values, the .integrand 

is nearly singular at some point. The integrand starts at zero, goes 

through a strong positive peak and through a strong negative peak,, 

Essentially, one has to subtract the effects of the two large and nearly 

2 

equal peaks. As Y or X gets closer to 1, and as the loss factor Z 
decreases, the peaks get narrower and steeper until the integrand resembles 
the doublet function. 

Now one can appreciate the numerical difficulty. The standard 

technique of singularity subtraction is not applicable because of the 

complicated form of the integrand. One cannot even find analytically the 
zero crossing of the integrand. An integration scheme with moderate-? . 
sized uniform steps might only randomly sample the nearly-singular region 
or overstep it completely. If the steps were small enough to accurately 
sample the peaks, there would be an excessive number over regions where . 
the function is relatively flat. Therefore, we have chosen a variable-step 
Romberg integration scheme [9]. 

The Romberg quadrature scheme is based on Richardson extrapolation 
of successive trapezoidal-rule estimates of the area under the integrand. . 
Suppose one needs the integral of some function f(x) over the interval 
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REAL Part, z 
X = .1, Y 2 = 
Z = .05 



REAL Part, z = 0 
X = .3, = .5, 

Z = .05 


REAL Part, z = 0 
X = . 7, Y 2 = .5, 
Z = .05 



REAL Part, z = 0 
X = .9, Y2 = .9, 
Z = .05 




IMAG Part, z = n 
X = .3, Y 2 = .7, 
Z = .05 


! r\ 



REAL Part, z = 0 
X = 1.5, Y 2 = .5, 
Z = .05 



IMAG Part, z = 0 
X = 1.9, Y 2 = 1.9, 
Z = .05 


Figure 4. 


Normalized Integrand in the Numerical Integration of 

E ,, for Various Plasma Parameters and Distances 
zf 

(k Q = .576, H = tt/2, n = 21) 
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[a, b] . Let h = a - b. The trapezoidal-rule estimate of the area for 
k 

2 equally spaced segments is 


T „ k ‘ pr' [ 1 £(x o ) + £(lt i> + ‘ ' + f(x 2 k_ 1 ) + 1 £ < lt 2 k )1 


(4.1) 


where 


x n ' a+ <n - £) 7 ' 


Extrapolations T ^ can be generated by the formula 


T "t = — - — (*“ T i , - T . ) . 

nk ^n _ v n-l,k+l n-l,k 


(4.2) 


These results can be tabulated as a triangular array 


00 




01 

T 

10 



02 

T 

11 

T 

20 


03 

T 

12 

T 

21 

T 

30 


(4.3) 


The elements of the first column are the trapezoidal-rule estimates; 
each other element is generated from the element to the left and the 

I . 

element above the element to the left. The accuracy of the estimates 

becomes better as one moves down and to the right. The procedure stops 

when the last two elements of a row are sufficiently close. Recall .that 

3 

the trapezoidal rule has an error term on the order of h . The application 

3 

of (4.2) to Tqq and Tq^ subtracts the h error term leaving only an . 

error on the order of h^ for T^. Substituting . (4.1) into (4.2) shows 

that is exactly the Simpson rule estimate; similarly, T£q is. the 

til 

third Newton-Cotes formula. One might suspect that T^g is the k 
Newton-Cotes formula, but fortunately that is not true. High- 
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order Newton-Cotes formulas suffer severely due to round-off error 
while Romberg’s method converges very well. 

Due to the particular integrand in this problem, we have 
partially followed modifications suggested by Miller and Burke [10], 

This modification is applicable to any numerical integration where 
the function is well-behaved except in a small subinterval. The original 
interval is divided into m subintervals and Romberg's method is applied 
to each. Convergence is checked by means of a test ratio T^. If 



< T 

r 


(4.4) 


only the midpoint of the subinterval needs to be computed and T^g can 
be used as the area of the subinterval. If (4.4) does not hold, two 
additional points in the sub interval are computed and the program checks 
if 



T . 

r 


If (4.5) is true, T^g is used as the estimate of the area. However, if. 
(4.5) is not satisfied, the program subdivides the interval and applies 
Romberg's method on that new sub interval father than subdividing all 
intervals as the unmodified version would require. Thus, no estimate 
higher than T^q is computed for an interval. This scheme is efficient 
since the midpoint of the new subinterval has already been calculated. . 
If one had tried using an adaptive Gaussian quadrature scheme, one could 
not have used any previous calculations. A distinctive feature of our 
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method is that our program saves all calculations made* Thus when 

a subinterval converges, one already has the function values for the end 

points and midpoint of the next sub interval to be integrated. To be . 

practical, one must limit the maximum number of subdivisions allowed 

and the minimum size of a subinterval. 

The exact error of this method is hard to estimate precisely. 

One may calculate an error estimate as the sum of the differences between 

the two best area estimates for each subinterval. This is a very 

pessimistic guess of the error. Even so, in most cases the relative 

error calculated in this way will be far less than T although there is 

no assurance that this is always true. The usual value used for T 
-4 

is 10 . 

A very interesting pathological case exists for this modified 
integration scheme. The T^g estimate term of (4.3) is exact for a 

g 

polynomial of degree up to 5. When the program tried integrating y= x 

over [0, 1] as a test problem, it kept subdividing the first interval 

until it had reached the maximum number of subdivisions allowed. The 

problem is analogous to integrating a parabola with the trapezoidal rule. 

2 

Consider approximating the function y = x with a straight line over 

[0, e] . The true value of the integral over the interval will be 

13 13 

— e but the integral of the approximation will be y e . Clearly the 

absolute error of the approximation decreases as e decreases but the 

relative error does not. When a maximum number of subdivisions is specified, 

there is no practical problem. Absolute error on the first interval is 

small enough, and for any interval not containing the 0, the relative error 

converges. In practice, this pathology occurs when the function behaves as 

a high-order polynomial. 
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The Romberg method has been described for the case of a real 
function. For a complex function, such as in this problem, we have 
required that the test ratio be satisfied for both the real and 
imaginary parts of the estimates before going to the next interval. 

The integration of is done slightly differently than formally 

indicated. The integrands of an ^ ^ n 2 are a ^ded together before 

integrating. This not only eliminates the need for two separate 

applications of Romberg's method, but also simplifies the handling 

of terms common to P , and P 

nl n2 
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CHAPTER V. STABILITY OF METHOD OF MOMENTS 

After many trial runs, we wished to test the stability of our 
method. Even if the models used were perfect and electric fields could 
be precisely analytically formulated, there would still be computational 
errors that would.be transmitted to the final answers. Knowing the 
sensitivity of the final answers to both matrix manipulation errors and 
errors in computing the fields allows one to estimate how accurately 
fields must be calculated to get results that are reliable. Early 
results obtained from field calculations with deliberate perturbations 
showed a noticeable sensitivity. The stability theory developed here 
applies not only to this particular plasma problem, but also to any 
problem using the method of moments. 

This sensitivity analysis will use norms of complex vectors and 
matrices to estimate error. Because of the wide variety in definitions 
of matrix norms, we wish to develop our analysis in detail and carefully 
define our matrix and vector notation for complex matrices and vectors. 

Our definition will follow Forsythe and Moler [11]. 

. T T 

All vectors x will be assumed to be column vectors. A and x 

will indicate the transpose of matrix A and vector x, that is, the 

H 

original matrix or vector with rows and columns interchanged. A and 
x will indicate the matrix or vector that is the complex conjugate of 
the transpose of the original matrix or vector. A matrix that is 
equal to its own transpose conjugate is said to be Hermitian. 

Next we define vector norms. The Euclidean norm of a complex 
vector x is defined as 

(5.1) 
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where x^, • • *»X n are components of x. For real vectors this 

norm is the usual concept of the length of a vector. We can also 
define vector norms 


n 


= l i*-l 

i=l 1 

(5.2) 

= max | x . | . 

(5.3) 

l<i<n 1 



These norms are not quite as geometrically obvious as the Euclidean norm. 
Figure 5 shows all possible unit-length vectors in two dimensions for 
each of the three norms. Mathematically, one can define an infinite 
number of norms, but these three are the most useful. All vector 
norms must satisfy the following three properties. 

| | cx | | = | c | • | | x | | for any real scalar c , (5.4) 

| |x| | 2 0 and | |x| | = 0 if and only if x = null vector, 

(5.5) 

||x+y|| < ||x|| + | | y | | (triangle inequality) . (5.6) 

Next we define a matrix norm for square matrices A as 

I I Ax | | 

f | A | | = max — - — - (5.7) 

x^O | | x M 

or equivalently 

| | A | | = max | | Ax | | . (5.8) 

11 * 11=1 

Each of the three vector norms thus produces a matrix norm that is said . 
to- be -compatible with the generating vector norm. Because_of this 
compatibility, these matrix norms satisfy the three vector norm 
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properties (5.4), (5.5), and (5.6), when x and y are replaced by matrices. 
From the definition (5.7) one can see 


| Ax | | < | | A | | • | |x | | for all A and x. (5.9) 


The main advantage of the two non-Euclidean norms is their ease of 

computation. It can be shown [12] that 

n 


I I A I I x = ma * l 

l<j.<n i=l 



I |A| | m = max l 
l<i<n j=l 



(5.10) 


(5.11) 


Finally we can introduce a scalar function of a matrix that 
measures the degree of ill-conditioning for certain operations with 
this matrix. 

Define the condition number of square matrix A as 


cond (A) = j | A | | • | |A 1 1 | . 


(5.12) 


For real matrices, the condition number with respect to the Euclidean 
norm can be calculated 


cond 2 (A) = /l 1 /X n (5.13) 

where A^ and A^ are the largest and smallest eigenvalues in absolute 

T 

magnitude of the symmetric matrix AA . If matrix A is symmetric 


cond (A) 



where X, and X 
1 n 


are the largest and smallest eigenvalues of A. 


' (5.14) 
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If matrix A is complex, (5.13) can be applied where one uses the 

H 

largest and smallest eigenvalues of AA which will be Hermitian and 
thus have real eigenvalues. 

For the case of real symmetric matrices and the Euclidean norm, 
one can most easily understand the application of the condition number. ; 

One can think of premultiplication by a matrix as a vector function in 
which components of the vector in the direction of the eigenvectors of 
the matrix are stretched or shrunk by the respective eigenvalues. 

Consider y = Ax as a mapping of vectors x onto vectors y. It is 
geometrically clear that y = Ax maps the n-dimensional sphere of all 
unit length x's onto an n-dimensional ellipsoid with semi-major 
axis A. and semi-minor axis A . One can see: that |A, I/Ia I is the 
ratio of possible distortions a vector x can undergo. Also note that 
the condition number shows the loss in accuracy for computation of components 
of x that are not in the direction of the eigenvector associated with the 
largest eigenvalue. The geometrical interpretation with other norms is 
not as clear. 

It is to be stressed that cond (A) is a far better measure of ill- 
conditioning than the smallness of the determinant A. Even if A is 
normalized, the smallness of det(A) has a weak relation to the degree 
of ill-conditioning. Supporting examples are given by Forsythe and 
Moler. 

Now we will indicate how the condition number can be used to 
measure sensitivity of a solution to a linear system. Consider the 
problem 

Ax = y (5.15) 

with unknown vector x. Suppose now that y is precisely known but that 
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A and, therefore, x are subject to uncertainty. Then 

(A + 6A) (x + 6x) = y (5.16) 

where 6A and Sx are uncertainties in A and x. Forsythe and Moler show 
from Equations (5.4) to (5.9) that 

J I fix | | _ I I «SA | | 

< cond (A) . (5.17) 

I |x + fix | | | | A | | 

Thus we can relate relative uncertainty in A to relative uncertainty 
in x by a function of A itself. It should be noted that (5.17) is the 
sharpest inequality possible; i.e., the equality will occur for some. .. 

<5A. Recall that any of the three norms can be used in (5.17). Occasionally 

• ||6A||. (5.18) 

To make practical calculations of the condition number, the 
condition number based on the infinity norm has been used. Equation 
(5.10) shows that the norm of A is the sum of the magnitudes of the 
elements of the row that gives the maximum sum. For the special 
form of the impedance matrix Z, this is always the center., row. For 
the inverse matrix Z \ in general, one must check all the rows above 
and including the center row. However, if k Q is small and the current 
distribution is essentially triangular, then the center row will have . 
the maximum sum of the magnitude of the elements. Multiplying the .norm 

of Z by~ the _ norm‘ of Z ^ gives the condition number which determines 

maximum error for the current on any segment. It should be noted 


it is useful to write (5.17) as 


6x 


: iia’ 1 ! 


x + 6x 
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that since Z. and Z ^ are symmetric, cond (Z) = cond 1 (Z) . Thus, if 
one estimates error in Z as | | 6Z | |^, one can use the same condition 
number to measure relative error in current I as | | 61 1 | / | 1 1 + 61 | | ^. 

Although the use of the Euclidean norm is geometrically clearer, 

it requires far more computations. First, one must form the Hermit ian 

matrix ZZ before finding eigenvalues. This multiplication, in itself, 

3 

takes on the order of n multiplications. Rather than finding .all the 
eigenvalues, one needs to find only the largest and smallest. A 
computer program has been written to find the largest eigenvalue by 
iteration. The iteration scheme operates on the principle that 
successive premultiplications of a vector by a matrix orients the 
resulting vector in the direction of the eigenvector associated with the 
largest eigenvalue in magnitude. Each iteration takes on the order of 
n operations and the convergence of this scheme is so poor that the 
convergence can be hastened by a linear extrapolation of two successive 
vectors. Finding the minimum eigenvalue is even harder. All of the 
eigenvalues of a matrix can be shifted by subtracting the maximum 
eigenvalue from the diagonal, thus making the former minimum eigenvalue 
the new maximum eigenvalue in magnitude. This, however, does not 
appear to be an efficient numerical procedure because small eigenvalues 
tend to be close together causing slow convergence. Even worse, more . 
accuracy is required so that when the largest eigenvalue is added to the 
final answer to obtain the true minimum eigenvalue, all significant 
figures are not lost. One can invert the matrix and find the maximum 
eigenvalue of that matrix. This eigenvalue will be the reciprocal of 
the minimum eigenvalue of the original matrix, but the inversion itself 
will take, on the order of n operations. All considered, using the 
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infinity norm, which is easy to calculate and gives the error in a 
Chebyshev sense, appears to be more practical. 
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CHAPTER VI. RESULTS 

Current distributions and antenna impedances have been calculated 

for various plasma parameters. Tables 2 through 9 and Figures 7 through - 

14 give the results for eight different plasma conditions. A graph 
2 

with axes X and Y is very useful for mapping plasma conditions. (See 

Figure 6.) The lines plotted on this graph correspond to various 

resonances in a lossless plasma. When dispersion surfaces are drawn 

over this plot, this display is known as a CMA diagram [13]. From this display 

one can see that the nature of the characteristic waves in the plasma 

changes abruptly as one crosses any one of these lines; hence, these 

lines form regions with similar properties. Each region has been 

sampled in the results presented. For our formulation, it is necessary 

to have some collision loss; therefore, for all cases, except for free 

2 

space which corresponds to X = 0 and Y = 0, Z has been arbitrarily 
set to .05. In each case H = 1, a = 1/250, and results are presented 
for three values of k Q . Since there is symmetry around the center 
feed, the graphs of current distribution show the current from the end 
point to the feed. The condition number based on the infinite norm, 
as described in the last chapter, is presented with each calculation. 

Balmain has developed an analytical expression for the input 
impedance using a quasi-static approximation which assumes a linear 
current distribution [14]. The impedance predicted by Balmain's 
formulation is presented for comparison. Where our method has determined 
a linear current distribution, the two values should agree. 

The first case presented is the degenerate plasma condition of 
free space. (See Table 2 and Figure 7.) One can see that the data 









63 


TABLE 2 . 

■ FREE-SPACE RESULTS (X = 0, Y 2 =0, Z = 0) 
(See Figure 7) 


e = 1.0 


e* = 0 


e = 1.0 
z 


IMPEDANCE (OHM)- 
COMPUTED 

ADMITTANCE (MILLIMHO) -r 
COMPUTED 

COND 

00 




.1600 - j 4491. 

7.937x10 + j .2227 

27.1 


-2 


10.39 - j 471.4 

4.671x10 +j 2.120 

35.8 

68.86 + j 16.52 

13.73 - j 3.295 

139 
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000 
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000 
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000 
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100 
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- 0 . 

750 

KQ- 

-1 . 

500 


MID POINTS OF SEGMENTS 


MID POINTS OE SEGMENTS 
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Figure 7. Current Distribution for Free Space (X=0,Y = 0, Z = 0) 






UNDERDENSE ELLIPTIC PLASMA RESULTS (X = .1, Y = .1, Z = .05) 
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2 

Figure 8. Current Distribution for Underdense Elliptic Plasma (X = .1, Y = . 
Z = .05) 


1 



UNDERDENSE ELLIPTIC PLASMA RESULTS (X = .4, Y = .5, Z = .05) 
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Figure 9. Current Distribution for Underdense Elliptic Plasma (X = .4, 
Y = .5, Z = .05) 




UNDERDENSE HYPERBOLIC PLASMA RESULTS (X 
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Figure 10. Current Distribution for Underdense Hyperbolic Plasma (X = .9, 
Y 2 = .9, Z = .05) 


UNDERDENSE ELLIPTIC PLASMA RESULTS (X=.5, Y = 1.5, Z= .05) 
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Figure 12. Current Distribution for Overdense Elliptic Plasma (X = 1.5, 
Y 2 = .5, Z = .05) 


OVERDENSE ELLIPTIC PLASMA RESULTS (X = 2.5, Y = .5, Z = .05) 
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Figure 13. Current Distribution for Overdense Elliptic Plasma (X = 2.5, 
Y 2 = .5, Z = .05) 


OVERDENSE HYPERBOLIC PLASMA RESULTS (X = 2., Y = 2., Z = .05) 
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Figure 14. Current Distribution for Overdense Hyperbolic Plasma (X = 2., 
Y 2 = 2 . , Z = .05) 
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agree generally with well-known results. The condition number for the 
impedance matrix Z is alarmingly high. For k Q = 1.5, which corresponds' 
to a dipole slightly longer than a half wavelength, the condition _ 
number is 140. This means that an error of .1 per cent in formulating 
the self term could cause "a 14 per cent error in the input impedance. 

One should be aware that for free space the. Green's function is much 
simpler and allows one to use better expansion functions in the method 
of moments. Mautz [15] and Klock [16] have done computer work for 
more arbitrarily shaped, thin-wire structures that agrees better, with 
measurements than our results for free space. However, condition numbers 
for some of these results are still high. 

Results for anisotropic plasmas are given in Tables 3 through 9 
and Figures 8 through 14. The validity of results varies with the 
region. Although there is not a unique wavelength in an anisotropic 
plasma, in several cases the current distribution looks like a 
frequency-scaled, free-space current distribution. Results appear 
best in regions’ which make an antenna look shorter , and worse in _ . 
regions which make an antenna appear longer, just as in free space 

i 

the stability was best for short antennas. For example, when H= 1 

and k Q =1.5, an antenna in free space is approximately half-wave 

2 5 

resonant, but an antenna in a plasma with X = 2 and Y = 2 is 
approximately full-wave resonant. When the current distribution is no 
longer linear, Balmain's approximation is no longer valid and there is 
no standard analytical approximation. 

Several regions require special comment. In the overdense hyper- 
bolic region, our results predict a negative radiation resistance for 
k Q = 7.5 and 1.5, which is clearly incorrect. The overdense elliptic 
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i 2 

region (X ! > 1, Y < 1) is more believable but has some unusual 

properties. The current decreases rapidly from the feed, making the 

current distribution curve concave upward, the reverse of the usual 

curvature of a sinusoidal distribution. The condition numbers indicate 

a much better conditioned problem than free space. Unlike other regions, 

the condition number decreases and the impedance agrees better with 

Balmain's approximation as k^ increases, although the current distribution 

is moving away from the assumed linear distribution. In general, regions 

separated only by -a hybrid frequency resonance show no drastic changes. 

There is little experimental work on current distribution which 

can be compared to our theoretical results. Some work was done by 

Snyder and Mittra [17], but they were unable to accurately determine 

the plasma parameters. An attempt was made to compare our results to 

work done by Ishizone et al. [18]. Their results did not indicate their 

antenna radius or a plasma loss parameter Z and so our computations 

were made with an arbitrary H/a ratio of 250 and a Z parameter of .05. 

Their measurements are all for overdense plasmas (X from 4. to 50) but 

2 

in both elliptic and hyperbolic regions (Y from 0 to 1.04). In all cases 

their current distributions have standing-wave characteristics, often 

with the current increasing at the feed. However, for all of their 

cases for which we have made computations, our current distribution 

simply decreases rapidly from the feed similar to our computed results 
2 

for X =2.5, Y = .5. Thus, there is rather poor qualitative agreement 
between our results and their measurements. 

In trying to understand why blatantly incorrect answers are 
s ome times obtained, one can note certain patterns between the 
condition number, the validity of the impedance, and the relative importance 
of the four terms of the Green's function E^. (See Tables 10, 11, and 12.) 
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As expected, the impedance matrix Z is best conditioned when it is 

most nearly diagonal, or when the singular terms which decrease rapidly 

with distance are dominant over the finite term. For the overdense 

elliptic plasma which gives good results as compared to Balmain's 

approximation, the field contribution on a segment distance 2H from 

the source due to the finite term is two orders of magnitude less 

than the total field, while for the overdense hyperbolic case X = 2, 

2 

Y =2, which gives negative radiation resistances, the finite term 

is an order of magnitude larger than the total field. The hypothesis 

that the difficulties encountered in some regions are due to an 

incorrect evaluation of the term, either through invalid approximations 

or poor numerical quadrature, is strengthened by the tendency of results - 

to get worse as k Q increases. In the singular terms, k Q appears only 

in the coefficient, while in the finite term, k is a factor in the 

o 

neglected terms of a Taylor series. A method is currently being 
investigated to evaluate fields for which the finite term is significant 
by an asymptotic expansion of the transformed Green's function. 

Calculations were made to test the sensitivity theory discussed 
in the last chapter. Consistent with the high condition numbers, it 
was found that a 1 per cent perturbation of the diagonal term in the 
Z matrix can cause a 30 per cent change of the impedance both for 
free space and the hyperbolic overdense plasma. One should note that . 
we are unable to obtain the condition number of an error-free matrix. 

If errors have occurred which make the matrix more diagonal, a 
deceptively low condition will be calculated. In cases where large . . 
-error-must- be-present— such-as— f or the overdense-hyperbolic-plasma,_one 



cannot, claim that the results approximate the true results within 
definite limits when only the condition number of the matrix with 


errors is known. 
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CHAPTER VII. CONCLUSIONS 

Using the method of moments we have been able to calculate the current 
distribution on an antenna oriented parallel to the static magnetic 
field in a plasma. For most plasma regions, reasonable impedances 
are obtained that agree with the Balmain approximation when k Q is 
small. Calculations on the overdense hyperbolic plasma, however, 
give unrealizable results. The exact cause of this discrepancy is 
not well understood; work is continuing on a different method of 
calculating the fields that will be more accurate and faster for 
segments far from the feed. 

A method of quantitatively analyzing the degree of ill-rconditioning 
has been discussed. This method applies not only to this problem but 
to any application of the method of moments. In particular, it is 
suggested that the stability of free-space, thin-wire problems be > 
evaluated. This method can also be applied to inversion problems 
arising from remote sensing studies. 



APPENDIX (LISTING OF PROGRAM USED) 
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C PREPARES E* S TO BE INVERTED 

C LIMITS CHANGES 

COMPLEX E > El t EZ, E ZS, EZFS1 , T ( 55 ) , EZ FS2 , EZ FN1 t T 1130 25 ) ,F Z Z , A D , 7 IMP, 
l E Z F N 2 » U t ER ( 2 1 1 »YH 2U i tHOLO » T E Z ( 2 1 ) 

DIMENSION PM( 30 ) ,PP ( 30) , ZZ ( 30 1 , SM ( 2 ) , PM1 ( 3 , 30 ) ,PP 1(3,30) 

REAL SCI ( 1 2 > /lC. ,8., 7.5,6. , 5. ,4. , 3, , 2. 50,2.', 1.5^, 1. ,0. / 

REAL K0/.576/,PER(3)/0. ,,01,.l/ 

REAL SSS(4>/4. 53, 5*35,13.47,50./ 

R^AL SP(2I/-3. , 1. 50/, Kl(3)/« 1,0.75 , 1. 5/ 

P 1=3.141593 

H = l. 

A= 1 . / 25C « 

CALL CCPIPH.5,1. , -3 ) P 

C M IS NUMBER OF SEGNENTS 

B=SQRT(2. )*A _ _ 

M = 21 ' ‘ 

XM=M _ 

DO 1000 KJ0B=1 ,1 

READ (5,251 XC,YC2,ZC 

25 FORMAT (3F10.6) 

YC = SQ RT( Y C2 ) __ ___ 

SS=0. 

SS1 = 0. _ 

NPLOT = 3 “ ' 

C XC , YC , AND ZC ARE 3 PLASMA PARAMETERS (DON’T USE X, Y,Z) 

DL=H/(XM*2.) ‘ “ “ ” ' ~~ 

U = ( 1. ,0. )-(0. ,1. )*ZC 

E=( 1. ,6. )-U*XC/(U**2-YC**2> 

E1=YC*XC/(U* *2- YC**2 ) _ 

E 1 = E 1 / E 

EZ = ( 1. .0. I-XC/U 

EZ=EZ/E 

__Q0__11 I XX=1 , NPLOT _ 

K0=K1( IXXI 

PRINT 24 , K0 , U , XC , YC2 ,ZC, E , E 1 , E Z ♦ M 

24 FORMAT ( • 1 • ,'KO = • , F 6 . 3 , 5 X , • U= ", 2 G 1 6 .7, 5 X , • XC = • , G1 6.7, EX, • YC? = • , 

1 G16. 7,5 X, 1 ZC= * , 

1 G16.7/'0 , ,'E=',2G16.7,5X,'E1= , ,2G16.7,5X,*EZ=»,?G16.7,5X, , M=',I4) 
PRINT 4 

4 FORMAT ( 1H1 , *>1 ELD ELEMENTS'/* NO. • » T 1 7 , • EZ S • , T4 3 , * E Z F S 1/ F Z F S 2 ' V 

1 T75 ,*EZFN1 I ,T109, , SUM | ) 

C EC YL FOR EZS, E2DL1 FOR EZFS1 

KK = M I NO ( 30, M) 

DO 2 I=1,KK 

__X=I-1 _ 

Z=X*2. *0L 

CALL E2DL1 (KG,E,E1 ,EZ,B,DL, Z , F ZZ , E ZF S1 , E ZF S2 ) _____ 

CALL FCYLtKO, E , E 1 , EZ , A ,DL , Z , EZS I 

CALL EZFN(K0,E,E1,EZ,DL,Z, EZFN1 , ER ( I ) , Y1 ( I ) ) 

T< II=EZS*EZFS1*EZFN1*EZFS2 
2 WRITE (6,511 l,EZS,EZFSl ,EZFN1,T(I) ,EZFS2 

51 FORMAT (1H , 1 3 ,6 El 5. 4 , 2E 1 8. 7/ 1 H ,33X,2E15.4) 

I F ( M , L T . 3 I I GO TO 12 
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C E2DL1 ; 

KK=M I NO ( 50 * M I 

DO 5 1 = 3 1, KK 

X= I- 1 

Z=X*2._*DL ’ _ 

CALL E2DL1(K0,E,E1 ,EZ,B,DL,Z*EZS,EZFS1 ,E"ZFS2I . 

CALL EZFN(KO,E,E! ,EZ,DL,Z, EZ FM 1 , ER { I) , Y 1 ( II) 

T( I )=EZS*EZFS1+EZFN1+EZFS2 

5 WRITE J6 f 51 l__L*EZS_, EZFS1 _ ,EZFN1Vni) ,EZFS2 

12 WRITE (6,1041 

104 FORMAT ( 1H0 ,• INTEGRATION AND ERROR • //_*N0. • , T 1 5 , • UNNPRMAL I 7E0 EZFN! 
T' ,T56, ' ERROR* ) 

C0ND-0» ; . 

00 71 1=1, M 

WRITE J6.ll I,Yl(IJ,ER(n 

3 FORMAT (IH , I 3 ,5 X , 2 El 5. 7 , 5 X , 2F 1 5, 7 ) 

_D0 6 J = 1 ,M_ _ 

K=M *( I-1I+J 

KK= I AB S ( J~* I > *1 

6 Tl(K»=T(KK» 

MLrl ABS(M/2tl- I» +1 

C0ND = C0ND+CABS( T (Ml » > 

71 CONTINUE 

WRITE (6,721 CONP 

72 FOR MAT ( 1H C ,5X NORM BE FOP E INVER SION =» , E 1 Jb. 7 ) 

CH0LD=C0ND 

CALL LINEQ(M ,T1» 

1 HOLD =0 
COND I T=0 • 

DO 73 1=1, M 

COND = Q« 

DO 74 J= 1 , M 
KK=M*( I-! ) ♦ J 

74 COND=COND+CABS (T1 ( KK I ) 

IF (CONO .LE.CONDIT) GO TO 73 
I HOL 0= I 

COND I T =COND 

73 CONTINUE 

WRITE (6,75) C OND I T , I HOLD 

75 FORMAT (1H ,5X,'N0RM AFTER INVERSION = • , E 1 6. 7 , 5 X , • ROW • , I ? ) 
COND=COND*CHOLO 

WRITE (6,76) COND 

76 FORMAT (1H , 5 X , ' FST I MATE OF CONDITION =» ,E16.7) 

PRINT 7 

7 FORMAT C • 1 • I 
KK=M**2 -2 

9 FORM AT ( • ,14,31 5X,2E16, 7) ) 

PRINT 15 

15 FORMAT {*€*///* NOT MULT BY -1/ ( 2DL )**’ 2» ) • 

DO 14 1=1, KK, 3 

14 PRINT 9, I ,T1 ( T I ,T1 ( I + 1 ) ,T1 ( 1+2 > 

C 

M2=M**2 
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F=-1./I2.*DLI»«2 

DO 13 1=1, M2 
13 TUI )=F*T1( I I 

MM=< ( M2-1 >/2+l > ' : 

AD=TU MM l*1000« 

ZIMP=1./T1(MM) 

Ml=((M-lt/2> * M 

M3=(M+l)/2 

M4=Ml*l 

M5 = M 1 ♦M3 

DO 17 I=M4,M5 
XX=CABS(Tl <n ) 

X Y=AT AN2 ( A I MAG ( T 1 ( I ) >,REAL(Ti( III) 

17 T 1 C I » =CMPLX< XX, XY ) 

PRINT 21 

21 FORMAT ( * O' ///• ', 10X , 'CURRENTS “ TMAG "AND PHASE*' ) 

DO 20 1=1, M3 

20 PRINT 22 »T1< I + Ml I 

22 FORMAT ( ' ' , IPX , 2G1 6« 7 ) 

PRINT 23,AD,ZIMP 

23 FORM AT ( ' O' /// ' • , 10X , ' ADM I TTANCE= *,2G16.7,' MILL I MHOS ' / 

"l • ' * 10X, • IMPEDANCE 3 * , 2GT6TT , • OHMS'l “ 

C START PLOTTING 

100 M4=M~3 + 1 

D O 30 1 = 1, M3 

PM (1*11= REALITK I ♦M'll I 
IF (PM( 1+1 >.GT,SS> SS = PM(IU! 

PM1 1 I X X , 1 + 1 I =— PM ( I +11 
PP(J+1I= A IM A G ( T 1 ( I + Ml I ) 

30 PPK IXX,I + l>=-PP(i + U 

PM ( 1 ) =Q« 

PM1( IXX,1*=0.0 

P PIll =c , 

PP1 ( IXX, 1 )=0, 

11 CO NTINUE _ . 

SS=SS/4* 

DO 53 K= 1 , 7 

IF ( SS.GE. 1. > "GO TO 54 

53 SS^SS*iO. 

54 K = K-1 

P£ 55 1=2^12 ___ •_ _ 

IF (SS.GT.SCK U I GO TO 56 

55 CONTINUE , 

56 SM(2)=SC1( I-ll/10*7*KU.E-6 

SM(1)=Q. 

C START PLOTTING 

CAL L CCP1 PA _ 

CALL CCP1PL (4.5, 2. ,-3 ) 

CALL CCP1PL (0. ,0« ,3 ) 

CALL CCP1PU0. ,6. ,21 

C_ALL FACTOR j^5 »_ 

CALL CCP2SY (-7*9, .5,, 3, * H= ',90.,3> 

CAJ.L C CP3NRJ 0 ^_ L 0 . , - ?_3 , _H ,9 0 ) 








CALL CCP2SYI-7.4,, 5 ,. 3, *X= «,90.,3 
CALL CCP3NRI0. ,0. ,-.3,XC,90. ,3) 

CALL CCP2SY (-7. 0 , . 5 , . 3 , « Y 2 =«,90.,5) 

CALL CCP3NRI0. ,0. 3,YC2,90. ,3> 

CALL CCP2SYI-6.6,.5,.3 , * 1 = »,90.,3) 

CALL CCP3NRI0. ,0. 3,ZC,90.,3) 

YP = . 5 

XP=-6. 2 

DO 52 K= 1 , N PLOT . 

L=K 

I F IK. EQ» 3 ) L=C __ 

CALL CCP2SY(XP,YP»«3» , K0=' ,90. ,3) 

CALL CCP3NR ( 0. .0. 3 , KI ( K > , 90. , 3 ) 

CALL CCP2SY ( 0. ,0. 3 , • . . . • , 90. , 3 > 

CALL CCP2SY ( 0. , 0. , 1 5 , L , 90. 1) 

52 XP=XP+.4 

CALL FAC T0R(1. ) 

CALL CCP5AXI0. ,0. ,*MAG OF CURRENT • ,14,4.000,180. , SM ) 

CALL CCP2SY I .3,1. 370,. 15, ' MID POINTS OF SEGMENTS « ,90. , 22 ) 
CALL CCP2SY ( .3 , 5. 800 ,. 1 5 , • FEED • , 90. , A \ 


XXI=I-1 

ZZ < I >=12./IXM+1. )»X XI 

CALL CCP2SY <0. , ZZI l ) ,.07,13,90. ,-l > 
DO 46 K=1 , NPLOT 






91 


SUBROUTINE LINEQ(LL,C) 

C 2/5/71 C KLEIN CHANGED ARITH IPS TD LOGICAL IFS 

COMPLEX C(l) ,STOR,STO,ST,S 
DIMENSION LRI58) 

DO 20 1 = 1, LL 
LR( I 1 = 1 

2 0 CONTINUE 

Ml =0 

DO 18 M= 1 , LL 

K=M ' ' 

DO 2 I=MjJ.L 

Kl=Ml+i 

K2=Ml+K 

IF ( <CARS<C<K1 M-CABS(C(K2H > .LE, 0. » GO TO ? 

6 K =1 

2 CONTINUE 

LS = LR(M> _______ 

L R ( M I = L R ( K J 

LRI K > = LS 

K2=Ml+K 

ST 0R=C(K2> 

J1=0 

DO 7 J=l,LL 
K1'=J1 + K 

K2=Jl+M 

STO=C ( Kl ) 

CIKl >=C(K2) 

C(K2 )=STO/STOR 

Jl=Jl + LL _ _ _ 

7 CONTINUE 

K 1 =M1 ♦ M 

C( Kl I = 1. /STOR 

D0_1_1 1 = 1, LL __ 

I F ( I — M I 12,11,12 

12 K l = M l ♦ I 

ST = C < K 1 } 

C ( Kl ) =0« 

J1=0 

DO 10 J = 1 , LL _ 

Kl = J 1 ♦ I 

K2 = J1»M 

C( K1»=CIK1 l-C(K2)*ST 

J1 = J 1 -*-LL 

10 CONTINUE 

___ 11 CONTINU E 

Ml =M 1 +LL 

13 CONTINUE 

J1 = 0 

DO 9 J = 1 , L L 

IFU-LR(J) .EG. 

14 LRJ=LR( J> 

J2=( LRJ-1 >*LL 

_ _21 DO 13 1= 1, LL 


0) GO TO 8 






K2= J 2* I 

K1=J1+I 

S=C( K2 ) 

C(K2»=C(K1J 

C ( K1 ) =S 

13 CONTINUE 

LR ( J > = LR ( LR J t 

LR(LRJ»=LRJ 

IF(J-LR(J) .NS* 0) GO TO 14 
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SUBROUTINE EZEN(KO , E t El , E Z , DL , Z , E IF N1 , ERROR , Y 1 ) 

REAL KO 

REAL»8 XY/1DO/ 

COMPLEX E, El, EZ»X,EZEN1, ERROR, Yl 

REAL* 8 P12/1. 570796326 8/ 

COMPLEX *16 Y,0APHI1 , DAP H 10 , DHPHI 2 , ERR 

X=DHPH 10 (DBLE(KO)tDCMPLX(DBLE(REAL(F ) ) t DBLE ( A IMAG( E >> > » 

1 DCMPLX(DBLE(REAL(E1> ) , DBLE ( A I MAG( F 1 ) ) > , 

2 OCMPLX ( DBLE (REA L ( £ Z ) ) , DBL E( A IMAG( E Z >)> , 

3DBLE ( DL > t DBLE ( Z ) I 

CALL ROMBR G_( ODO_,_P 1 2 , D APHI1 V B,Y ,1 0- 4, 10-3 , ERR > 

ERROR=ERR 

Y1=V 

K=OR EALZ ( DHPH I? ( XY ) I 

x = y . ; 1_ 

EZFN1=(0* ,15, )*KG**2*X 

RETURN ' ' 

END 
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SUBROUTINE E2DL1 ( KO , E, E 1 1 E Z ♦ A , PL ♦ Z t EZ S . E ZFS1 ,E 7FS2) 

C THIS SUBROUTINE CALCULATES FIRST THREE E TER M$ USING 2 DL MODEL 

COMPLEX E,E1 ,EZ,EZS,EZFS1 ,EZFS?,A1?,ZAP, Z AM, ZAPH, ZAMH,ZF?,ZF? 

C RETURNS COMPLEX NUMBERS 
C THIS REPLACES SUBROUTINE E2DL __ 

REAL KO 

A12=£Z*A*»2 

ZPDL = Z +DL 
ZMDL=Z- DL 
ZPOL2=ZPDL**2 
ZMDL2=ZMDL**2 
ZAP = ZPDL**2 *A12 

Z AM= 2 MDL **2 +A12 

ZAPH = CSQRT( ZAP I " ~ ' ~ " 

ZAMH-CSQRT(ZAM) _____ _ 

Z F 2 = CLOG! ( ZPDL + Z APHI/ ( ZMDL+ZAMH) I 
ZF3=ZPOL /ZAPH - ZMDL/ZAMH 
C EZS SIMPLIFIED 2/7/70 

E Z S = ( 0 . ,15. )/ (K0*DL*E >*( ZPDL/ ( ZAP* ZAPH I -ZMDL /( Z AM* ZAMH) > 
EZFS1=(0. ,-7. 5l*<KO/DL>*( (2. ,0. I*ZF2-ZF3) 

IF ( REALIEZI .NE. 1. .OR. AIMAGIEZI _.NF. .OR. RFAL (Ell 

1 .NE. 0. .OR. A I MAG ( El I .NE. 0.1 GO TO 1 
EZFS2= ( 0. ,0. ) 

GO TO 2 

1 li^ALOG ( ( ZPO L+S QRT ( ZPQ L2 +A** 2 I > / ( Z MD L ♦ SORT ( ZMDL? »A* *2 ) ) I 

E Z F S 2 =Y 0 . « 1 5. I *E"l**2 / ( E Z- ( 1. , 0. I > •* ?* ( KO/DL I ' ( T1 + ZF2* ( ■ 5 ) ; ' 

1 ( -2. ,0. I + ZF3M .5,0. I’M ( 1. »0. I-EZII 
C 1 (1.- EZI I 

2 RETURN 
END 
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SUBROUTINE ECY L ( KO , E , El t E Z 1 A , DL , Z ♦ E ZS ) 

C THIS EVALUATES EZS USING ELLIPTICAL INTEGRALS 

CO MPLEX E»El > EZ»EZSfX,A 12» KP2 > KM2 

' REAL KO 

C THIS FUNCTION STATEMENT EVALUATES ELLIPTICAL INTEGRAL 
ELm = ((((A4*X + A3)*X + A2 >*X + A1 l*X + l. ) 

1 - ( ( ( ( B4* X+B3 ) *X *B2I*X +B1 )»X)»CLOG(X) 

DATA Al, A2.A3, A4/ . 4432 5 1 4 , 6. 2 6060 IE- 2, 4, 757384E-?,!. 736506E-?/ 
„ DATA B1,B2 , S3 t B 4/« 24 9983 7 1 9 . 200 1 80E-2 , 4.T 69698E- 2 , 5 , 2 6449 6E-?/ 
ZPDL=Z+DL 

ZMDL = Z-DL __ 

A12=EZ*A**2 

KP2 = 4« *A12 / ( ZPDL**2+4« * A12 ) 

KM2=4» ♦A12/(ZMDL**2+4.*A12> 

EZS= (0» t !5« )/( (K0»3. 141 593«DL*A )*£* C SO RT ( EZ \ _ 

1 (C SORT (KP2I * E L ( ( 1 .76. ) -KP2I/Z PDL- 

2 C SQRT <KM2 )*EL( ( l.tO. ) -KM2 ) / Z MD L > 

RETURN 

END 


96 


SUBROUTINE ROMBR G ( A , B. TTT , NM I N , V , TR , HM I N , ERR I 

C CONVERTED TO COMPLEX 

IMPLICIT CQMPLEX*16IA-G,0-W,V-ZI • 

IMPLICIT REAL*8IH,XI 

R EAL *8 TR , A, B ,DREALZ,Q1MAGZ 

C TNT MUST BE DECLARED EXTERNAL IN CALLING PROGRAM 

C ADAPTIVE INTERVAL ROMBERG INTEGRATION ROUTINE 

C NM IN IS MIN NUMBER OF SUBINTERVALS 

DIMENSION FF2 1 10 I . FF4 I 10 ) , HF I 10 I 

NSUB=10 

KODE =0 _ __________ 

HM AX = ( 8-A) /NMIN 

ERRg ( ODO.OOOI 

Y = ( 000 »CDO I 

X4 = A 

F4=TNT(A» 

DO_l_J_Jrl,NMIN __ : ; _ 

C INITIALIZE LARGEST INTERVAL POSSIBLE 

X0 = X4 

F0 = F4 ~ ' 

X4=X0+HMAX 

F4=TNT ( X4 I 

X2 =X 0 + * 5^H MAX 

LEVEL=0 

F2 = TNT ( X2 ) 

H = HM A X 

C START LOOP _ __ 

4 W = FC + F4 

C H. FAC TOR ED OUT OF TIN, Ml TERMS 

TOQ= W/2D0 

WW=W+2DQ*F2 

TO 1 = WW/ 4D0 

T10= ( 4D0*T0l-T00 I /3D0 _ 

IF ( DREALZ ( T 1 0 I .EQ.O.ODOI GO TO 2 0 

IFIDMAX1 (DABS! (OREALZITIO I -DRE ALZ I TCI > I /DREALZ I Tl^ I ) , 

I DABSKDIMAGZITlOI-DIMAGZITOin/DIMAGZITioi I I . i_E. 

2 TRI GO TO 8 

GO TO 21 

20 I F ( DABS.(_(_DIMAGZ IT 101-01 MAGZI TOi I ) /DIMAGZI T10I I . LE. TR I GO TO R 

21 CONTINUE 

F 1 = T N T I X 0_t_« 2 5 * H I _ . 

F3=TNT(X4-.25*H) 

TC2 = ( WW+2DQ*(F1+F3 I 1/803 

T11=(4DO*T02-T01 I/3D0 

T 2 0 = 1 1 6 DC * T 11 - T 1 Q J. / 1 5 DO . _ 

IF I DREALZ I T20 I .EQ.O.ODOI GO TO 22 

IFID MAX1 I DABS! I DRE AL Z I T20 ) -DREAL Z I T 1 11 I /DR EALZ I T2" 1 > I ♦ 

l 0ABS((DIMAGZ(T2°I-DIMAGZIT11I l/DIMAGZITll) I I 

2 > l e, „ tri go to . 2 : 

GO TO 23 

. 22-1 F I _D ABSI I DIMAGZIT20 i-DIMAGZ l Tllll/DI M I AG Z I _T_1 1 U L E ,.TR )__GCL TO 2 

23 CONTINUE 

I F I LEVEL . LT. _NSUB_. AND. H .GT.HMI NI GO TO 3 
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K0D£=1 

GO TO 2 

C SUBDIVIDE 

3 H=H/2. 

LEVEL = LEVEL +1 
FF2I LEVEL) =F3 
HF( LEVEL )=H 
FF4< LEVEL I = F4 
F4=F2 
F 2 = F 1 
X4=XO+H 
GO TO 4 

8 Y=Y+TIO*H 
ERR=ERR+ (TIO-TOI )*H 
GO TO 9 

C UNSUBDIVIDE 

2 Y=Y+T20*H 

ERR=ERR-MT20-T1I »*H 

9 IF ( LEVEL .EQ. 01 GO TO 1 
F0=F4 

F2=FF2 ( LEVEL l 
F4=FF4 ( LEV EL ) 

H=HF( LEVEL) 

X0=X4 

X4=X0+H 

LEVEL=LEVEL-1 

G O TO 4 

1 CONTINUE 

IF ( KODE .EQ, 1 ) PRINT 5 

5 FORM AT ( ' ACCURACY SPECIFIED CANNOT "BE OBTAINED* ) 

RETURN 

END 



O O 0,0 o o 
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COMPLEX FUNCTION DHPHI0+16I KO, E, El , EZ ,DL , Z ) 

IMPLICIT REAL*8( A-H,0-Z > 

EXTENSIVE CHANGES 7-28-71 TWO INTEGRAND HALVES SUMMED AND COUNTER 
FOR TIMES CALLED ADDED 

CHANGES MADE 5/1/71 HC ELIMINATED AND INCLUDED IN HPHT 
CHANGES MADE 4720/71 ABOUT SQPT OF N1 
DON'T DO UNIAXIAL CASE 

CHANGES MADE 4/15/71 HPHI NO LONGER NEFDS HCEO 
COMP LEX* 16 Nl,Ul,U2,COEF,H,Ei,EZ,E,N12,ST,nSCP,N22,DUM,F12,A,p.,C, 

1 H 1 » Dl » T N T , D HP H I 2 

REAL*8 KC 

1=0 

DHPHIO=(CDO,ODO) 

RETURN 

ENTRY TNT( PHI I 
1=1*1 “ 

SP = DS I N ( PH I ) 

SP2=SP**2 
C P = DCOS ( PH I I 

I F ( CDABS(El) .EQ. COO GO TO 1 

C H FOR El -,= 0 _ 

E 1 2= El **2 “ 

CP2=CP**2 
A = SP2 + FZ* CP 2 

B = EZ»I 100»CP2>*(1D0 -E12 I»SP ? 

C=( 100-E12 »*EZ 
DSCR=CDSQRT( B**2-4D0*A*C I 
N12 = (B-DSCR)/(2r)0*A> 

N22=( B+DSCP I / ( 2D0* A ) 

DUM= (ODO«-2DO)*SP*SP2/(A**2*(N12-N22) I 

D l = lDr.*E12*CP? 

H=OUM* ( 1 DO — E 1 2-N 1 2 * D 1 ) 

HI = — DU M* ( 1 DO — E 12— N22*D1 I 
GO TO 2 

C H FOR JE1»0_ _ 

1 N12=(1C0,CD0) 

N22=N12 

H = (0D0,2DD |*SP*SP2 
H.1 = ( ODOjDDO ) 

C CALCULATING S IN < >/ ( ) TERM 

2 S T =C DSQRT(N12*E) 

IF (DIMAGZISTI « GT, 0001 ST=-ST 

Ul = ( ST*KC*CP ) 

U2=U1*DL 

IF J _ CDABSIU2 1.LT. 1.D-3I GO TO 

COEF=CDS INIU2 > /U2 
GO TO 4 _ 

3 COEF=( 1D0,0D" » 

4 TNT = H*CDEXP<-( O DQ I DO )» U1»DAPS( Z ) I *COEF + ST 

ST=CDSQRT(N22*E) 

IFIDIMAGZIST > »GT. GDC I ST = -ST 
Ul=( ST*KO*CPf 
U2=U1*DL 


IF ( C0ABSIU2 ).LT. l.D-3) GO TO 6 

COEF=CDSIN<U2)/U2 

GO TO 8 ; . 

6 COEF=< lDO,ODO) 

8 TNT=H1* CDEXP(- <0 00» 100)»Ul»DABSm )*COEF* S T ♦ TNT 

9 RETURN 

ENTRY DHPH 1 2 ( X ) 

DHPH i2=DCMPLX( OBLE ( FLOAT { I » ) ,000) 

RETURN 


END 
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